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Abstract 

We compute analytically the statistics of the Renyi and von Neumann 
entropies (standard measures of entanglement), for a random pure state 
in a large bipartite quantum system. The full probability distribution is 
computed by first mapping the problem to a random matrix model and 
then using a Coulomb gas method. We identify three different regimes 
in the entropy distribution, which correspond to two phase transitions in 
the associated Coulomb gas. The two critical points correspond to sudden 
changes in the shape of the Coulomb charge density: the appearance of 
an integrable singularity at the origin for the first critical point, and the 
detachement of the rightmost charge (largest eigenvalue) from the sea 
of the other charges at the second critical point. Analytical results are 
verified by Monte Carlo numerical simulations. A short account of some of 
these results appeared recently in Phys. Rev. Lett. 104, 110501 (2010). 



1 Introduction 

Entanglement plays a crucial role in quantum information and computation as a 
measure of nonclassical correlations between parts of a quantum system [T] . The 
strength of those quantum correlations is significant in highly entangled states, 
which are involved and exploited in powerful communication and computational 
tasks that are not possible classically. Random pure states are of special interest 
as their average entropy is close to its possible maximum value [5J [3] . Taking a 
quantum state at random also corresponds to assuming minimal prior knowledge 
about the system [4]. Random states can thus be seen as "typical states" to 
which an arbitrary time-evolving quantum state may be compared. In addition, 
random states are useful in the context of quantum chaotic or nonintegrable 
systems [SJ [HI E]- 

There exist several measures for quantifying entanglement [5]. For a bi- 
partite quantum system, the entropy (either the von Neumann or the Renyi 
entropies) is a well-known measure of entanglement. For a multipartite system, 
the full distribution of bipartite entanglement between two parts of the system 
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has been proposed as a measure of multipartite entanglement [S] ■ The distribu- 
tion of entropy in a bipartite system is thus generally useful for characterizing 
entanglement properties of a random pure state. 

Statistical properties of observables such as the von Neumann entropy, con- 
currence, purity or the minimum eigenvalue for random pure states have been 
studied extensively [3 [3 HOI HH H3 H3 HH [T31 HQ H3- In particular, the av- 
erage von Neumann entropy is known to be close to its maximal value (for a 
large system). In contrast, few studies have addressed the full distribution of 
the entropy: only the distribution of the purity for very small systems [13] and 
partial information on the Laplace transform of the purity distribution for large 
systems [TU] have previously appeared in the literature. 

Our purpose here is to compute the full distribution of the Renyi entropies for 
a random pure state in a large bipartite quantum system. In particular, we show 
that the common idea that a random pure state is maximally entangled is not 
quite correct: while the average entropy is indeed close to its maximal value |H 
[3J , the probability of an almost maximally entangled state is in fact vanishingly 
small. This statement requires to compute the full probability distribution of 
the entropy, namely its large deviation tails, which is one of the goals achieved 
in our paper. 

The calculation of the Renyi entropies' distribution proceeds by mapping 
the entanglement problem to an equivalent random matrix model, which de- 
scribes the statistical properties of the reduced density matrix of a subsystem. 
We can then use Coulomb gas methods borrowed from random matrix theory 
We identify three regimes in the distribution of the entropy, as a direct conse- 
quence of two phase transitions in the associated Coulomb gas problem. One 
of those transitions is akin to a Bose-Einstein condensation, with one charge of 
the Coulomb gas detaching from the sea of the other charges - or equivalently 
one eigenvalue of the reduced density matrix becoming much larger than the 
others. 

This paper is a detailed version of a short letter that was published re- 
cently 18 . It thus contains all explicit formulas for our results and details 
about analytical proofs and numerical simulations as well as new results, es- 
pecially for the third regime of the distribution (see below), the von Neumann 
entropy and the maximal eigenvalue of the density matrix. 

The plan of the paper is as follows. In section [2j we describe precisely our 
model of bipartite quantum system for the direct product %a <8> T~Lb of two 
Hilbert spaces Ha and Hb ■ In section [3J we analyze the distribution of the 
eigenvalues Aj of the reduced density matrices of the two subsystems. In partic- 
ular, we compute the average density of eigenvalues and explain the Coulomb 
gas method that we also use later for computing the distribution of the Renyi 
entropy S q — j^\nY, q where S g = y\ Xf. In section 01 we compute the full 
distribution of £ 9 for a large system. We find two phase transitions in the asso- 
ciated Coulomb gas, and thus three regimes for the distribution of S 9 . In section 
[3 using results from section [4j we derive the distribution of the Renyi entropy 
Sq as well as the distribution of the von Neumann entropy (case q — > 1) and the 
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distribution of the largest eigenvalue (q — > oo). Finally in section [BJ we present 
results obtained by Monte Carlo numerical simulations that we performed to 
test and verify our analytical predictions. 



2 Random bipartite state 

In this section, we set the problem of bipartite entanglement for a random pure 
state. We first describe a bipartite quantum system, introduce then measures 
of entanglement (the von Neumann and Renyi entropies) and give finally the 
precise definition of random pure states. 

2.1 Entanglement in a bipartite quantum system 

Let us consider a bipartite quantum system A® B composed of two subsystems 
A and B of respective dimensions N and M. The system is described by the 
product Hilbert space Wab = 'Ha^'Hb with N = dim (Ha) and M = dim (Kb)- 
Here, we shall be interested in the limit where N and M are large and c = ^ 
is fixed. We shall take N < M, i.e. c < 1, so that A and B play the role of the 
subsystem of interest and of the environment, respectively. 

Let \ip) be a pure state of the full system. Its density matrix p = \tp){^}\ 
is a positive semi-definite Hermitian matrix normalized as Tr p — (iplip) = 1. 
The density matrix can thus be diagonalized, its eigenvalues are non-negative 
and their sum is unity. Subsystem A is described by its reduced density matrix 
PA = Ti'b [p] = J2a B =i( aB \p\ aB )> wnere \ aB ) is an orthonormal basis of Hb- 
Similarly, B is described by ps = Tta [p\- It is easy to show that the reduced 
matrices pa and ps share the same set of non-negative eigenvalues {Ai, Xn} 



Any pure state can be written as — E a =i Xi , a \^ A ) ® \ aB ) where 

\i A ) ® \a B ) is a fixed orthonormal basis of Wab- The singular value decomposi- 
tion of the matrix x^ a permits to recast the previous expression in the so-called 
Schmidt decomposition form: 



where \mf) and \pf) represent the eigenvectors of pa and pb, respectively, 
associated with the same eigenvalue A^. 

The representation ([1]), namely the Schmidt number n$ of strictly positive 
eigenvalues, is very useful for characterizing the entanglement between subsys- 
tems A and B. For example, let us consider two limiting cases: 

(i) If only one of the eigenvalues, say A^, is non zero then Aj = 1, ng = 1 
and the state of the full system \ip) — \mf) ® \pf ) is a product state, which is 
said to be separable. The system is unentangled. 



with Y,i=i A» = 1. 



N 




(1) 
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(ii) If all the eigenvalues are equal (Aj = 1/N for all j), ns — N and is 
a superposition of all product states. The system is maximally entangled. 

A standard measure of entanglement between two subsystems A and B 
is the von Neumann entropy of either subsystem: SVn = — Tr [pa In Pa] — 
— Xi ln\i, which reaches its minimum when the system is unentanglcd 
(situation (i) above) and its maximum IniV when the system is maximally en- 
tangled (situation (ii)). Another useful measure of entanglement is the Rcnyi 
entropy of order q (for q > 0): 

(2) 

which also reaches its minimal value in situation (i) and its maximal value 
In TV in situation (ii). As one varies the parameter q, the Renyi entropy inter- 
polates between the von Neumann entropy (q — » 1 + ) and — lnA max (q — » oo) 
where A max is the largest eigenvalue of the reduced density matrices. 

2.2 Random pure states 

A pure state is called random when it is sampled according to the uniform Haar 
measure, which is unitarily invariant. Specifically, a random pure state is defined 
as \ip) = J2iLi Elli x i,a V A ) ® \ aB )' where \i A ) <£> \a B ) is a fixed orthonormal 
basis of Hab and where the variables {xi^ a } are uniformly distributed among 
the sets of {xi, a } satisfying the constraint J^i a \ x i,a\ 2 = 1 (normalization of 
\tp)). Equivalently, the probability density function (pdf) of the N x M matrix 
X with entries Xi yCC can be written 

P(X) cx 6 (Tr(AAt) - l) oc e -^ Tr ( xxt ) 8(Tr(XX*) - 1) , (3) 

with the second equality showing that the pdf can also be seen as a Gaussian 
supplemented by the unit-trace constraint. 

In the basis \i A ) of *Ha, the reduced density matrix of subsystem A is simply 
given by pa = XX^ . In general, when X is a N x M Gaussian random matrix, 
i.e. P(X) cx e^2 Tl '( xxl ) (iid Gaussian entries x^ a that are real for a Dyson 
index j3 = 1, complex for (3 — 2), the N x N matrix XX* is a Wishart matrix 
whose distribution of eigenvalues is [15] : 

N 

i— 1 i<j 

The Vandermonde determinant Yli<j — •M^ ma kes that the eigenvalues are 
strongly correlated and they physically tend to repel each other. 

The major difference between the matrix pa = XX* in the quantum problem 
and a standard Wishart matrix stems from the unit trace constraint Tr [pa] — 1- 



N 

.1=1 
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The constraint is to be included in the distribution of the eigenvalues of pA, 
which is given [3j [11] by: 



N 

P(\i,...,\n) = B MiN 5(^2^-1) J]Af (M - JV+1) - 1 JIlAi-A^, (5) 

i i—1 i<j 

with (3 = 2 (the x^ a are complex) and the normalization constant Bm,n com- 
puted using Selberg's integrals [TT] : 

B = T(MNf3/2)T(l + l3/2) N 

M ' N n?Jo l n(M-j)l3/2)r(l + (N-j)P/2)' 

The presence of a fixed trace constraint (as in Eq. (JSJ) is known to have impor- 
tant consequences on the spectral properties of a matrix J5UJ HI] • We will see 
that in the present context also, the fixed trace constraint does play an impor- 
tant and crucial role. In particular, this constraint is directly responsible for a 
Bose-Einstein type condensation transition that will be discussed in the context 
of the probability distribution of the entanglement entropy. 

Since the eigenvalues Ai of pa are random variables for a random pure state, 
any observable is a random variable as well. Statistical properties of observables, 
namely of various measures of entanglement such as the von Neumann entropy 
[51155]. G-concurrence [T5J, purity [TUIHB] or minimum eigenvalue [TH [TCI [T01 [T7] . 
have been studied extensively. In particular, Page [3] computed the average von 
Neumann entropy in the limit M > N 3> 1: (Svn) ~ In TV — M*. He also 
conjectured its value for finite TV and M (which was proved later [55]). In 
contrast, there have been few studies on the full distribution of the entropy, 
except for the purity S2 = Af whose distribution is known exactly for small 
TV (2, 3 and 4) [15) . For large TV, the Laplace transform of the purity distribution 
(generating function of the cumulants) was studied recently [10] for positive 
values of the Laplace variable. However, when inverted, the previous quantity 
provides only partial information about the purity distribution. 

Here, we compute analytically the full distribution of the Renyi entropy S q 
(defined in Eq. (J2j) or equivalently of = YliLi K = CX P [(1 — l)S q ], for large 
TV. As for the von Neumann entropy, the average value of the Renyi entropies is 
close to their maximal value In TV (maximal entanglement) : (S q ) ~ In TV — z(q), 
where z(q) > (for q > 0) is independent of TV for large TV. For example, for 
M w TV and q = 2, we have z(q = 2) = In 2. However, we show below that the 
probability that S q approaches its maximal value In TV is again very small. 



3 Distribution of the eigenvalues of pa 

The eigenvalues of the reduced density matrix pa are distributed according to 
the law in Eq. ([5]). Given this joint distribution, the first natural object to study 

is the average spectral density pn.mW = j? X)i=i — ^*))- This average 
density Pn,m W dX also gives the probability to find an eigenvalue between A 
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and X + dX (the one-point marginal of the joint distribution) . For finite (N,M), 
this average density was computed first for (3 = 2 [23l [24] and very recently 
for f3 = 1 25 . However, these formulae involve rather complicated special 
functions and taking the asymptotic large N, large M limit is nontrivial. Here 
we will take a complementary route which is well suited to derive exactly the 
asymptotic limit. We will take the limit N — > oo, M oo but keeping their 
ratio < c = N/M < 1 fixed. For the spectral density, we will henceforth use 
a shorthand notation pat(A) = Pn,n/cW- We will show that for large N the 
limiting form of Pat (A) can be obtained easily via using a Coulomb gas approach. 

Due to the unit trace constraint J2i=i ^« = 1; the typical amplitude of the 
eigenvalues is X typ ~ for large N. Since X typ ~ 4= (and pat is normalized to 
unity), we expect (as will be proved below) that the average density for large 
TV has a scaling form: 

Pn (A) « N p* (XN) . (7) 

Using the Coulomb gas method explained in subsection l3.1[ we find an exact 
expression for the rescaled density p*{x): 



P* ( x ) ~ TTTir V x - Li L 2 - x , (8) 

and 



Ittcx 

2 - i— x 2 



where the right and left edges read L 2 — c (y^ + ; -^l = c ~ l) 
we recall that c = iV/M < 1. 

For c — 1 (N M), Li = 0, L2 = 4 and the rescaled density reduces to: 



In Fig. [I] plots of the rescaled density p* (x) and comparisons to the shape of 
the rescaled density for a standard Wishart matrix are shown for c = 1 and 
c= 1/3. 



3.1 Computation of the rescaled density: Coulomb gas 
method 

The goal of this section is to prove Eqs. ([7]) and §8$ for the average density of 
states. The joint distribution of the eigenvalues in Eq. (JSJ can be interepreted 
as a Boltzmann weight at inverse temperature /3 

P(Ai, X N ) cx exp {—/3E [{A 4 }]} , (10) 

where the effective energy is given by 

N 

E[{Xi}} = -'yJ2lnXi-^2ln\X i -X j \ with ^ A, = 1 . (11) 

i— 1 i<j i 

Here, 7 = M ~^ +1 _ _j_ ^ jyCL_£) f or l ar g e AT. The logarithmic binary inter- 
actions correspond to the Coulomb repulsion in 2 dimensions. The eigenvalues 
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p*(x) 




X 

Figure 1: The rescaled average density p*(x) of the eigenvalues for the density 
matrix of a quantum subsystem. The rescaled density is defined by Pat (A) m 
N p* (\N) for large N (see Eq. ©) and is plotted for c = § = 1 (red solid 
line) and c = 1/3 (blue dashed line). The density is compared with the rescaled 
average density of Wishart eigenvalues (random matrix theory) : p^y(x) defined 

by Pn W « jr P*w (w) ( see ^ C3) P lottcd for c = f = 1 ( red solid line ) 
and c — 1/3 (black dotted line). The different dependencies on c for p*(x) and 
Pw( x ) ma ke that, even after their different rescaling in N, the two distributions 
are equal only for c = 1. 
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can thus be seen as charges of a 2D Coulomb gas, repelling each other electro- 
statically. The charges are confined in the segment 1 > A,; > for all i and they 
are also subject to an external logarithmic potential (with amplitude 7). 

The mapping from random matrix eigenvalues to a Coulomb gas problem 
is well-known in random matrix theory and has been recently used in a variety 
of contexts that include the distribution of the extreme eigenvalues of Gaus- 
sian and Wishart matrices [26l [27j [28l [29] , purity partition function in bipartite 
systems [TU], nonintersecting Brownian interfaces [30 , quantum transport in 
chaotic cavities |31j . information and communication systems |32) . and the in- 
dex distribution for Gaussian random fields [33 , 34 and Gaussian matrices [35] . 
Here, we use similar methods yet the problem is quite different due to the con- 
straint K — 1- First, the scaling with TV (for large N) differs from standard 
Wishart matrices. Indeed, X ty p ~ l/N in our problem of entanglement whereas 
AjL ~ N for a Wishart matrix. However, the effect of the constraint J^i = 1 
is not just the rescaling of standard Wishart results by a factor of l/N 2 as it 
may seem. It turns out that the constraint has more serious consequences and 
leads to fundamentally different and new behavior (including a condensation 
transition which is absent in Wishart matrices) that we will demonstrate. 

Configurations of the eigenvalues are characterized by the density p(A, N) = 
N^ 1 Yli=i $ — ^0- F° r large N, the eigenvalues are expected to be close 
to each other and their typical amplitude is X typ ~ jj- We introduce then a 
rescaled variable x ~ 0(1) as x — XN . The corresponding density is p(x) = 

N" 1 Eili 5 - W), so that P( X , N ) = N p(XN) = N p(x). 

The effective energy in Eq. (jlip becomes in the continuous limit (large N) 
a functional of the density p. To the leading order in N, the effective energy 
reads E [{AJ] = N 2 E [p] + 0(N), where 

— — \ J dx p{x) In x — — J J dxdx' p(x)p(x') In \x — x'\ 

+ ^0^ dx p(x) - l\ + fjti (^J dxxp(x)-lj. (12) 

The Lagrange multipliers po and p\ enforce respectively the constraints J p = 1 
(normalization) and J dxx p(x) = 1 (unit trace) . 

The joint distribution of the eigenvalues is given by the Boltzmann weight 
P(Ai, Xn) oc exp{-/3N 2 E[p] + 0{N)} for large N. This distribution is 
highly peaked around its most probable value p* which is thus also the mean 
value of p: p*(x) = iV" 1 £i=i<<5 (x - X.N)). Hence, the average density of 

states is the continuous density p* that minimizes the effective energy: = 

9 p=p* 

0. From Eq. (fl2|) we get the saddle point equation for p*: 

fl-c\ 

dx' p*(x) In \x — a; I = ^0 + Mi x — — ^ — hix . (13) 

V 2c J 
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Differentiating with respect to x leads to the integral equation: 



dx 



Pi 



2c 



(14) 



where V denotes the principal value. This singular integral equation can be 
solved by using a theorem due to Tricomi [35] that states that if the solution 
p* has a finite support [L±,L2], then the finite Hilbert transform defined by the 
equation F(x) = V J^ 2 dx' p JfJ can be inverted as 



p*(x) 



■ny/x - LavLi 



c-v 



dx' \Jx' — L\\fL^ 



F(x') 



, (15) 



where the constant C fixes the integral of p* via f L ^ dx p*(x) = C. 

In Eq. (fT4| . F(x) = fi\— (^r) -■ Physically, the average density is expected 
to be smooth and thus to vanish at L\ and Li (bounds of its support): p* (L\) = 
= p*{L-2). These two constraints fix the value of L\ and L>2- The other two 
constraints J p* = 1 and / xp* = 1 give the value of the constant C in Eq. (|15|) 
and the Lagrange multiplier p\ in Eq. (1141) . Finally, inserting the expression of 
p* in Eq. (JT3J) for a special value of x (say x — L2) gives po. Imposing all these 
constraints, we finally get: 



P*( x ) = o y/x- Lx \/L 2 



2ncx 



(16) 



with Li.a = c (l T \J\) (where c = N/M). We also find C = J p* = 1, 

pi = 1/ (2c) and no = - (^) +2 (l - ±) In [1 + ycj + Finally, the average 
density in the original variable A is given by pjv (A) = N p* (A7V), where p*{x) 
is given in Eq. P^|) . 



3.2 Comparison with Wishart eigenvalues 

For Wishart matrices, it is known that the average density of the eigenvalues is 
given, for large N and fixed c = N/M, by the Marcenko-Pastur law [37] : 

P%W*jf with p* w (x) = J-J X -LY sjLj - x , (17) 

with the right and left edges given by LY = (l + y^j an d ^Y = ^1 — yc) ' 

As expected, the scaling with N is different: A]^ p ~ N for a Wishart eigen- 
value, whereas the unit trace constraint makes that Xt yp ~ 1/N for an eigenvalue 
of the quantum density matrix pa- 

For c = 1, the two edges L}^ = 0, iJf = 4 and p* w {x) = p*(x). However, for 
a general c < 1 the rescaled densities are not quite the same (even though they 
have the same shape): P\y(x) — cp*(xc). Figure [T] shows a comparative plot of 
Py/{x) and p*(x) for c = 1 and c = 1/3. 
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4 Distribution of S g = ^ A* for q > 1 and c = 1 



This section is somewhat long as it contains the bulk of the details of our 
calculations. Hence it is useful to start with a summary of the main results 
obtained in subsections 4.1-4.3 as well as the main picture that emerges out of 
these calculations. Readers not interested in details can skip the subsections 
4.1-4.3 and get the main picture just from this summary. 

In this section, we compute the full distribution of H q — '^2 i X q , and thus 
of the Renyi entropy S q = In (S g ) /(l — q) for large N. We take for simplicity 
M « N, i.e. c = 1, but our method can be easily extended to c < 1 as well. 
For simplicity, we will also restrict ourselves to the case q > 1. However, our 
method is also easily extendable to the case < q < 1. Since = 1 

and x — > x q is convex for q > 1, we have iV 1-9 < S 9 < 1 (or cquivalently 
In TV > S q > 0). The lower bound S 9 = TV 1 "" 3 corresponds to the maximally 
entangled case (situation (ii) in subsection 12. 1|) . when Xj = 1/N for all j: the 
entropy is S q = InN. The upper bound E g = 1 corresponds to the unentangled 
case (situation (i) in subsection I2.1|) when only one of the Xi is non zero (and 
thus equal to one): the entropy is zero. 

The scaling Xt yp ~ 1/N implies that E 9 ~ N^ q for large N. Let s = 
T, q iV 9_1 be the rescaled variable s ~ 0(1). In figure [2l a typical plot of the 
probability density function (pdf) P (S 9 = iV 1-9 s) is shown: the distribution 
has a Gaussian peak (centered on the mean value s = s(q)) flanked on both 
sides by non-Gaussian tails. We show below that there are two critical values 
s = Si(q) and s — S2(q) separating three regimes I (1 < s < si(<?)), II (si(q) < 
s < s 2 {q)) and III (s 2 (q) < s). 

At the first critical point Si(g), the distribution has a weak singularity (dis- 
continuity of the third derivative). At the second critical point S2(q), a Bose- 
Einstein type condensation transition occurs and the distribution changes shape 
abruptly (first derivative is discontinuous in the limit N —¥ +oo). These changes 
are a direct consequence of two phase transitions in the associated Coulomb 
gas problem, more precisely in the shape of the optimal charge density. The 
schematic plot of the distribution of E g (for large N) in Fig. [2] clearly shows 
the three regimes I, II and III and the discontinuity of the derivative at ,s = s 2 
(transition between II and III). 

More precisely, the probability density function of S 9 for large N and q > 1 
displays three different regimes: 



The exact mathematical meaning of the "~" sign is a logarithmic equivalence : 




for 1 < s < si(g) ; 




for s > s 2 (q) ■ 



(18) 
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1 si(g) s 2 (q) s lsi(g) Js 2 (q) s 

Figure 2: Schematic distribution of S 9 = Yli = ^V 1-9 s as a function of s for 
(very) large iV. Panel (a) shows the shape of the pdf of E g , while (b) shows 
the shape of the rate function — lnP(E g = N 1 ^ q s). Two critical points si(q) 
and S2(q) separate three regimes I, II and III, characterized by the different 
optimal densities shown in figure [3] The maximally entangled state s = 1 is at 
the extreme- left of the distribution, well spaced from the mean value s(q). 



Pc(X,N) Pc (X,N) Pc (X,N) 




Li/N L 2 /N L/N ( t 

(a) (b) (c) 

Figure 3: Scheme of the optimal saddle point density p c of the eigenvalues (or, 
equivalently, of the Coulomb gas of charges) for (a) 1 < s < si(q) (regime I), (b) 
s i(q) < s < s 2(q) (regime II) and (c) s > S2(q) (regime III). In regime III, the 
maximal eigenvalue A max = t becomes much larger than the other eigenvalues, 
as shown by the isolated bump in (c). 
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(a) (b) (c) 

Figure 4: Scheme of the effective potential V(x) seen by the charges of the 
Coulomb gas (eigenvalues) for (a) 1 < s < si(q) (regime I), (b) si(q) < s < s 2 (q) 
(regime II) and (c) s > S2(q) (regime III). In regimes I and II, the charges are 
confined close to the minimum of the effective potential. In regime III, the 
potential is not anymore bounded from below. Therefore, one charge detaches 
from the sea of the other charges : the maximal eigenvalue becomes much larger 
than the other. 



In P( £ —N r ^ q s) 

/hv 2 — > $/(s) as N — > oo with fixed s <E [l,si(g)[ (resp. $77 for 

In P( £ —N r ^ q s) 

fixed s &]si(q), S2 (<?)[) and ffjvw/g — — ^ ^iii(s) as N — > oo with fixed 

s > s 2{q)- The rate functions $/, and (as well as si and s 2 ) are 

independent of N - but they depend on the parameter q. Explicit expressions 
of the functions and are given in Eqs. and (02]) for q = 2, and in 
Eq. (|47| for a general q > 1; an explicit expression of ^//j is given in Eq. (|50[) 
for a general g > 1 (and in Eq. (|5T|) for q = 2). As shown in figures [5] and [6] 
(resp. for N = 50 and ./V = 1000), we also did some Monte Carlo simulations 
(as explained in section [6|) and found that our analytical predictions agree very 
well with the numerical data. 

Regime II includes the mean value (S q ) rs A^ 1_g s(g), i.e. si(q) < s(q) < 
S2(q) for every q. The mean value is explicitely given by: 



<£,) « N^- S (q) With 5(g) = -^-^L ^. (19) 



r( g + i/2) 

y^T(q + 2] 

For large TV, the distribution of Y> q given in Eq. (fl8|) is highly peaked around 
its average (because of the factor TV 2 in regime II): the average value of S g 
coincides then with the most probable value, i.e. s{q) is the minimum of $//(s). 
The quadratic behaviour of $77 (s) around this minimum gives the Gaussian be- 
haviour of the distribution of S g around its average (and thus gives the variance 
of E g ). We get: 

P {Eg = N 1 ~ g s) P3 exp |-/3iV 2 ii_ for s close to s(q) . (20) 
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Therefore, the variance of E g is given by: 



2 



4 29 , ^V{q+l/2f 



VarE, = (Ej> - <£ 9 ) 2 » with o\ = —q{q if . (21) 

The distribution has a Gaussian peak flanked by non-Gaussian tails described 
by the rate functions $j (left tail) and ^in (right tail). Conversely, the rate 
function describes the middle part of the distribution, which includes the 
Gaussian behaviour in the neighbourhood of the average. 

In the limit N — > oo, si(q) and S2(q) do not depend on N and the second 
critical value S2(q) is actually equal to the mean value s(q) of s: 

T(q + S/2) fA(q + l)\ g r(g + l/2) ,„ ,„„, 

v/7rr(g + 2) V 3q J \AT(g + 2) 

However, for a large but finite N, 52(9, -/V) actually depends on N and is given 
in Eq. ([23]) below. 

The convergence in TV for the regimes I and II is very fast : the agreement 
between numerical simulations and analytical predictions in the limit N — > 00 
is very good already for N ~ 50. However, the second transition, between 
regime II and III, is affected by finite-size effects, that remain important even 
for N ~ O(10 3 ). Their main effect is a shift in the value of the critical point 
S2- The transition actually occurs at a value S2(q, N) that depends on N, is a 
bit larger than s(q) and tends slowly to s(q) as N — > 00. More precisely, the 
second transition occurs at s = 52(9, N) with 

y/qp(q-l)s(q) 

s 2 {q, N) « s(q) + J= .. fn _ lWf? J n for large but finite AT . (23) 



N (q-l)/(2q-l) 

For example, for q = 2, we have s(g = 2) = 2 and ,52(9 = 2, Af) 
2^1s^ f or i ar ge but finite N. 



N173 



The extreme left of the distribution corresponds to maximally entangled 
states: s — > 1 + means that A? = — » N 1 ~ q , that is the case where all the 
eigenvalues are equal and the state is maximally entangled (situation (ii)). As 
s — » 1, 3>/(s) tends to +00, thus the pdf P(T, q = N 1 ~ q s) tends rapidly towards 
zero. For example, for q = 2, we have P(S q = N 1 ~ q s) ps (s— 1)^ / 4 as s — >• 1 + . 
This implies that the probability of a maximally entangled configuration is very 
small (for large N). 

Similarly, the extreme right s — > +00 of the distribution corresponds to 
weakly entangled states. An unentangled state has indeed only one non-zero 
eigenvalue, A$, thus S = T, q = 1 (situation (i)). We can actually compute the 
expression of the pdf for the scaling E g = S with S ps 0(1) (S 3> s/N) and 

/ s/3N 2 /2 

< 5 < 1. For q = 2, we get: P (S 2 = 5) « M - x/SJ for AT -> 00 with 

S ps 0(1). For 5 — ?> 1 _ , the pdf of E 9 is again tending very rapidly towards 
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zero: unentangled states are highly unlikely. 

The three regimes in the distribution of £ q are actually a direct consequence 
of two phase transitions in the associated Coulomb gas problem, as we show 
in this section. We compute the probability density function P(E 9 = TV 1-9 s). 
The charges of the associated Coulomb gas see a different effective potential 
V(x) depending on the value of s, as shown by Fig. [4] 

• In regime I (1 < s < si), the potential V(x) has a minimum at a positive 
x and the charges accumulate near this minimum: the optimal density /0 C (A, N) 
describing the charges has a finite support over [Li/N, L2/N] and vanishes at 
Li/N and L 2 /N (see Fig.^a) andija)). 

• In regime II {s\ < s < s 2 ), the potential is minimum at x = 0, the 
charges accumulate close to the origin: the optimal density p c (X,N) describing 
the charges has a finite support over ]0, L/N], vanishes at L/N but diverges as 
l/\/A at the origin (see Fig.^b) andl[b)). 

• As s exceeds s 2 , the potential becomes unbounded from below; the right- 
most charge (maximal eigenvalue) suddenly jumps far from the other eigenval- 
ues: the charges are described in regime III by a density with finite support 
]0,C] and a single charge (maximal eigenvalue) well separated from the other 
charges: t ^> £ (see Fig. E^c) and|4](c)). 

4.1 Computation of the pdf of £ g : associated Coulomb 
gas 

In this subsection, we explain how we compute the pdf (probability density 
function) of T, q using a Coulomb gas method. The pdf of S 9 is by definition: 



The joint pdf of the eigenvalues P(Ai, A at) is given in Eq. and can be seen 
as a Boltzmann weight at inverse temperature /?, as in Eq. (1101) : 



where the energy E [{AJ] = -7 ^ j=1 In A t - Y^i<j ln ~ I ( with J2i ^ = 1) 



is the effective energy of a 2D Coulomb gas of charges. For large N, the effective 
energy is of order E ~ 0(N 2 ) (because of the logarithmic interaction potential). 
We can thus compute the multiple integral in Eq. (|2"4"]l via the method of steepest 
descent: for large N, the configuration of {A^} which dominates the integral is 
the one that minimizes the effective energy. 

For Eq. (f2"4")l we also have to take into account the constraint J^i ^1 = 
(delta function in Eq. (f24|0 . This will be done by adding in the effective energy 
a term f/ 2 (^ i A? — E g ) where fi' 2 plays the role of a Lagrange multiplier. Phys- 
ically, this corresponds to adding an external potential /j,' 2 A 9 for the charges. 





P(Ai, \n) oc exp {—(3E [{Xi}]} , 



(25) 
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For large N, the eigenvalues are expected to be close to each other and 
the saddle point will be highly peaked, i.e. the most probable value and the 
mean coincide. We will thus assume that we can label the Xi by a continuous 
average density of states p(X,N) — A^X^W^ ~~ ^)) = N p(x) with p(x) = 
N^ 1 J2i{fi( x ~ ^iN)) and x = XN . However, we will see that this assumption 
is not correct for large E g (large s): in the regime III, the maximal eigenvalue 
becomes much larger than the other eigenvalues. The maximal eigenvalue should 
then be treated on its own and be distinguished from the continuous average 
density. 

Let us begin with the case where the eigenvalues can be described by the 
density p[x). Then the pdf of S g can be written as: 

P(S 9 =N 1 ~«a,N) ex Jv[p] exp{~l3N 2 E s [p}} , (26) 



where the effective energy E s [p] is given by 
2 



E s [p] = — — / / dxdx' p(x)p(x') In \x — X 1 \ + po ( / dx p(x) — 1 



+pi ij dx x p(x) — 1 ) + P2 yj dx x q p(x) — sj . (27) 

The Lagrange multipliers po, p\ and pi enforce respectively the constraints 
J p = 1 (normalization of the density), Yli^i — 1 (unit trace) and = 
N 1 ~ q s (delta function in Eq.((24])). 

For large N, the method of steepest descent gives: 

P (E, = N 1 ^ s, N) cx exp {~(3N 2 E S [p c ] } , (28) 

where p c minimizes the energy (saddle point): 

5E ; 

T P 

The saddle point equation reads: 



= . (29) 

P=Pc 



/■OO 

/ dx p c (x') In \x — x'\ = po + p\x + p2% q = V(x) , (30) 
Jo 

with V(x) acting as an effective external potential. Differentiating with respect 
to x gives: 

V r dx' ^] =pi+q P2X q - x = V'(x) , (31) 



where V denotes the Cauchy principal value. The solution for a finite support 
density p c is given again by Tricomi formula as in Eq. (|15l) and yields the answer 
for the regimes I and II. 
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In these regimes, the pdf of E g is thus given by P (S g = TV 1 g s, N) ~ 
exp {— /37V 2 $(s)} where the rate function $(s) is equal to E s [p c ] up to an 
additive constant. More precisely, the normalized pdf reads: 

P P«- N S > N )~ fV[p]exp{-m 2 E[p}} ' (32) 

where E s [p] is given in Eq. (l27l) and E [p] is the effective energy associated to 
the joint distribution of the eigenvalues (without further constraint), as given 
in Eq. (|T2")) (we remind that c = 1 in the present section). The steepest descent 
for both the numerator and denominator gives: 

P ft - ., N) M ZWEW " - W'» ■ <33) 

with $(s) = E s [p c ] — E [p*] and where p* (resp. p c ) is the density that minimizes 
E [p] (resp. E s [p] ) . The density p* (x) is thus simply the rescaled average density 
of states given in Eq. ([9]) (for c — 1). Finally, we get 

$(s) = E s [p c ] - E [p*] = E s [ Pc ] - 1/4 . (34) 



4.2 Regime I and II 

Regimes I and II correspond to the case where the eigenvalues can be described 
by a continous density p(x), as explained above. In this case, we have seen that 
the pdf of Eg is given for large N by P (S, = iV 1 ^ s, N) rj exp {~/3N 2 <S>(s)}. 
In this section, we derive an explicit expression for $(s) = $>i(s) in regime I ie 
for 1 < s < si(q) (Eq. ([38]) in subsection 14. 2 .11 for q = 2) and $(s) = $//(s) in 
regime II ie for Si(g) < s < s 2 (q) (Eq. (|42|) for q = 2 and Eq. (|47| for a general 
q > 1 in subsection 14. 2. 2ft . 



4.2.1 Regime I 

The solution of Eq. (|3"Tj) is a density with finite support [Li, L2] where L\ > 0. 
As the density is expected to be smooth, we must have p c (L?) — and p c {L\) = 
at least for L\ > 0. As the eigenvalues \ are nonnegative, another possibility 
is that L\ = and p c (Li) ^ - this will be regime II. The first case, i.e. with 
L\ > and p c [L\) — 0, defines the regime I and is valid for 1 < s < Si(q) with 
si given in Eq. (|22|) . as we shall see shortly. 

In this subsection, we show that, for 1 < s < si(q) (regime I), p\ < and 
p.2 > 0, hence the effective potential V(x) defined in Eq. (|30[) has a minimum 

at a nonzero x: at a; = a;* = ^^77^ 91 > 0, as shown by Fig.QJa). The charges 
concentrate around this nonzero minimum. Thus the density of charges p c is 
expected to have a finite support over [L\, L2] with L\ > and to vanish at the 
bounds Li t 2 (see Fig. E2a)). 
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A simple case: q = 2 



Let us begin with the case q = 2, where we can find an explicit expression 
for the density p c and the pdf of the purity £ 2 = Yli ^1 = Tr [Pa] ■ 

We find the solution of Eq. pip for q = 2 by using Tricomi formula with 
F(x) — V'(x) (cf Eq. (IT51) ). The solution p c has a finite support [Li,L 2 ]. By 
imposing p c (L{) = = p c (£ 2 ) (regime I), we get: 

2/*: 



p c (x) = \fx - ~L\ y/ L 2 - x . (35) 



The optimal charge density is a semi-circle. At this point, there are six unkown 
parameters: the constant C in Tricomi's formula; the bounds of the density 
support Li and L 2 ; the Lagrange multipliers p$, p\ and p 2 . We also have some 
constraints to enforce. The two constraints p c {Li) = = p c {L,2), together with 
the three constraints J p c = 1, J xp c — 1 and J x 2 p c = s fix the value of the five 
parameters C, L±, L 2 , p\ and p 2 . We get po by inserting the final expression 
of p c in Eq. ([30]) for a special value of x, say x — L2. 

By imposing these constraints, we find C — J p c = 1, la, 2 = 1 T 2Vs — 1, 
^1 = ~ 2(3-1) ' ^ 2 = 4(/-i) and Mo = 5 In |s- 1| + 4(^1) - \- Therefore we have 



\/L 2 - x\/x - Li 

PC{X) = 27T( fl -l) ' (36) 

with Li ;2 = 1 =F 2\/s — 1. This solution is valid for L\ > 0, i.e. for s < 5/4. 
Thus, regime I corresponds to 1 < s < s\(2) with si(2) = 5/4. 

In this regime, we have p\ = — ^(s-x) ^ ^> M 2 = 4( s 1 _i) > 0> an d the effective 
potential V(x) = po + p\x + p 2 x 2 has a minimum for x — x* — 1 > 0. The 
charges concentrate around this minimum: they form a semi-disk centered at 
x* = 1 = (Li + L 2 )/2. The radius of the semi-disk R = 2\fs — 1 increases with 
s till L\ reaches its minimal possible value (for s = 5/4). 

Finally we compute the saddle point energy. Using the saddle point equation 
(Eq. J3Q])), we get E s \p c ) = -\ (p + pi + p 2 s) = — \ In (s - 1) + \, which gives 
the expression of $/(s) = E s [p c ] - E [p*] = E s [p c ] - \ (see Eq. ((34j))- The 
distribution of the purity T, 2 is thus given by: 

P(£ 2 = s/AT,A0cxexp{-/W 2 $ / ( S )} , (37) 
where the large deviation function is explicitly given by: 

<!>/(*) = -- In ( S -l)-~. (38) 



General case: q > 1 

The same qualitative behaviour holds for a general q > 1: in the regime I, 
the effective potential V{x) has a minimum at a nonzero x — x* > 0, the charges 
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accumulate around this minimum. The density p c has a finite support [L\,Li2\ 
with L\ > and p c (Li) = = p c (L 2 ). This regime is valid for 1 < s < si(q). 
The value of the critical point is determined from the analysis of regime II: we 
show that regime II is valid for s > si(q). Unfortunately, we were not able to 
obtain explicit expressions for p c and $/ in regime I for general q (the integral 
in the Tricomi formula for a general q seems hard to compute analytically). 

4.2.2 Regime II 

As s approaches s\(q) from below, the lower bound L\ of the density support 
tends to zero. As the eigenvalues are non-negative, L\ cannot be negative. 
Hence, regime I does not exist for s > s\(q). The critical value S\(q) is the 
onset of regime II, where the density p c has a finite support ]0, L] and vanishes 
only at the upper bound L (see Fig. |3jb)). We will see that regime II is valid 
for Si(q) < s < s 2 (g, N) where S2(q,N) is given in Eq. (j2"3"j) . 

Within regime II and for increasing s, p\ increases and becomes positive 
while pi remains positive. The effective potential V{x) = po + p\x + p2X q has 
thus a minimum at a smaller and smaller value x — x* that sticks to zero when 
pi becomes positive (see Fig.|4jb)). The charges concentrate close to the origin. 

A simple case: q = 2 

Let us begin with the simple case q = 2. We find the solution of Eq. (|3~Tj) for 
q = 2 by using again the Tricomi formula with F(x) = V'(x) (cf Eq. (fl"5|0. We 
are looking for a solution p c with finite support [0, L]. After imposing p c (L) = 0, 
we get: 

p c { x ) = lJhE^[ A + Bx] , (39) 

7T V X 

with A = pi + P2L and B = 2pi- 

There are five unkown parameters: the arbitrary constant C in Tricomi's 
formula; the upper bound of the density support L; the Lagrange multipliers 
Po, pi and p2- We also have constraints to enforce. The constraint p c {L) = 
together with the three constraints j p c = 1, j xp c = 1 and J x 2 p c = s fix the 
value of the four parameters C, L, p\ and pi- We get po by inserting the final 
expression of p c in Eq. (|30|) for a special value of x, say x = L. 

We find C = J Pc = 1, pi = 8(L - 3)/i 2 , p 2 = 4(4 - L)/L 3 and p a = 
ln(j) — \ — p\\. The upper bound of the support L is solution of the equation 
I? - 12L + 16s = 0. Hence L = 2(3 ± y/9 - 4s). Physically the density p c (x) 
must remain positive for < x < L. It is not difficult to see that this determines 
L: 

L = L{s) = 2(3 - V9 - 4s) (40) 

The upper bound L increases with s and matches smoothly regime I: L = 
2 = L2 at s = si(2) = 5/4. The solution of regime II, exists as long as 
s < 9/4. However, we shall see that there exists another solution for s > 2 that 
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Figure 5: Distribution of £2 = Yli '■ the figure shows the rate function $(s) — 

~ lnP ^pN 2 plotted against s for N = 50. Analytical predictions (red solid 
line) are compared with the results (blue points) of Monte Carlo numerical 
simulations (method 1, as explained in section [6]). Our analytical predictions 
consist of three regimes. For regimes I (1 < s < 5/4) and II (5/4 < s < 2), we 
have plotted the asymptotic expressions of the rate functions in the limit TV — >• 00 
given in Eqs. (|38l) and (|42l) . For regime III, we have plotted the analytical 
prediction for large but finite N, using for $ui(s, N) — $(iV, s/N) (see Eq. ([5"7|) ) 
the complete expression of E given in Eq. (|6ip and C and t (numerical) solutions 
of Eq. ([591 and Indeed, for iV = 50, finite- corrections to the asymptotic 
formula in Eq. (|5ip are important in regime III : the curve of the dominant 
behavior in N would not fit well the data and the complete expressions are 
needed. Note in particular that finite- N effects make that the transition between 
II and III is regularized and appears to be smooth. 
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is energetically more favorable. This latter solution will yield regime III. The 
solution of regime II is thus valid only for 5/4 < s < 2. 

We have seen that /ii = 8(L — 3)/L 2 and [i 2 = 4(4 — L)/L 3 . According to 
the respective sign of \i\ and /Z2, we distinguish three phases for the effective 
potential V(x) = Mo + /J,ix + ix^x 1 : 

• 2 < L < 3 (i.e. 5/4 < s < 27/16): Ml < and /i 2 > 0. The potential V(x) 
has a minimum at a positive x — x* = (— Ml) / (2 M2) = L{3 — L)/(4 — L) 
(as in regime I), x* decreases when L (or s) increases and reaches at 
L = 3 (see Fig. 4 (a)). 

• 3 < L < 4 (i.e. 27/16 < s < 2): /ii > and ^ 2 > 0. The potential 
is monotonic (increasing) on the real positive axis. It has an absolute 
minimum at x — (see Fig. 4 (b)). 

• L > 4 (i.e. 2 < s < 9/4): /^i > but ^ 2 < 0. The potential is not 
anymore bounded from below. It increases around the origin, reaches a 
maximum at x = x* = (fi±) / (—2 ^2) = L(L — 3)/(£ — 4) and decreases 
monotonically for x > x* to —00 (see Fig. 4 (c)). In this phase, the origin 
is a local minimum and the solution in Eq. f|39[) is metastable. There is 
actually a second solution in this phase, where one eigenvalue splits off the 
sea of the other eigenvalues. This second solution becomes energetically 
more favorable at s = S2 ~ 2 + jpja ■ The solution of regime II in Eq. (|39l) 
is thus valid only for s < 82- For s > S2, the second solution dominates: 
this is regime III. 

Finally, the distribution of the purity E2 in regime II is computed by the 
saddle point method: 

P (Sa = s/N, N) oc exp {-/3JV 2 $//(s)} , (41) 

where the large deviation function $77 = E a [p c ] — ~ = — i [/ii + /Z2S + Mo] — 4 
is explicitely given by: 

^ , n 1 , ^\ 6 5 7 

with L = 2 (3 — — 4s). For large N, this solution is valid for si(2) < s < 

s 2 (2,N) with si(2) = 5/4 and s 2 (2,N) w 2 + 2 as iV ->• +00 (as we 

shall see). 

At s = si = 5/4 (transition between regime I and II), the rate function 
$(s) has a weak nonanalyticity. It is continuous, $(5/4) = — | + and even 
twice differentiable: ^1 .,. = — 1 and 4-f I , = 4. However, the third 

ds ls=5/4 as ls=5/4 

derivative is discontinuous: €3\ = 4*^1 r/, = ~ 32 but C^l .... = 

ds- 3 \s=5/4~ as" l s =5/4 as" \s=5/4+ 

- , $ / J I _ .. = — 16. The minimum of $ is reached at s = 2 within regime II, 

ds d ls=5/4 ' 

which gives the mean value of the purity (£2) ~ 2/iV (as the distribution is 
highly peaked around its average for large N). 
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Figure [5] compares our analytical predictions for regimes I and II in Eq. 
and (|42p with numerical data (Monte Carlo simulations) : the agreement is very 
good already for TV = 50. 

General case q > 1 

We find the solution p c with finite support [0, L] of Eq. ((3T|) for q > 1 by using 
again the Tricomi formula with F(x) = V'(x) (cf Eq. ([TB]) '). After imposing 
p c (L) = 0, we get the expression of the density: 



AH jL-x , ^qLi- 1 r (q + \) jL-x „ f 1 1 3 1 a; 



^=tv— + ^ — r^rv— ^^.1-^,1-^,(43) 

where 2-F1 is a hypergeometric function 2^1(1, c, z) = Y^Lo ^fep^ 2 with 
(a) n = a(a+l)...(a+n — 1) denoting the raising factorial (Pochhammer symbol). 

Exactly as for q = 2, the constraints fix the unknown parameters. We obtain 
the Lagrange multipliers fj,\, p.2 and /iq as functions of L: 

8(1 + 9) 4q (1 + g) V^r(9) i-4 

Mi = 7 — ^ e4r 7— 2 — r and u 2 = ) ^ — 7 -r ■ (44) 

m (l-g)L 2 L(l-g) ^ (1 - <?) T{q + 1/2) Z^+i V J 

and /l*o = In (-j) + Mi 2g — ^ e u PP er bound L (which is a function of s) 
is given by the solution of the equation 

l + qj r(g+l/2) v 7 

For g = 2, we recover the simple expressions of the previous subsection. 

The function / : L -> L 9 + 4L?- 1 is increasing with L for < L < Lq 

with Lq = 4(1 + g)/g, and decreases for L > Lq. It is thus maximal at L = Lq, 
which implies that s cannot be larger than sq = s(L = Lq) in this regime. 
Hence, regime II is not valid for s > sq, where sq = So(q) = s(L = Lq) = 

r( 9 +l/2) ( 4(l+q) \ q 
2^T(q+2) \ q J ■ 

Moreover, it can be shown that, for L < Lq/3 and for L > Lq, the density 
p c (x) becomes negative for x close to the bounds (close to for L < Lq/3, close 
to L for L > Lq). This is not physical. Hence, L must belong to the interval 
[Lq/3, Lq]. Within this range, the function / is monotonic and it increases 
with L. It can thus be inverted and gives L as a single-valued function of 
s: L = L(s). This range [Lq/3,Lq] corresponds to si(q) < s < so(q), where 
si(q) = s(L = Lq/3) and s = s(L = Lq). 

Therefore regime II can exist only for s±(q) < s < SQ(q), where si(q) = 
(HiM) 9 and So{q) = (^)\~For , = 2, we recover 

Si(2) = 5/4 and sq(2) = 9/4. However, as in the q — 2 case, this regime is not 
valid anymore for s > s 2 (q, N) given in Eq. (|23|) . where a second solution starts 
to dominate (regime III). 
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Finally, we compute the pdf of E g as a function of L = L(s). We get the pdf 
by the saddle point method: 



P(E 9 = N^s^) oc exp{-/3iV 2 $ // (s)} , (46) 
where the large deviation function = E s [p c ] — j is explicitely given by: 

The function L = L(s) is the unique solution of Eq. (|45l) within the range 

S± < S < S2- 

Exactly as for q = 2, the parameter /i 2 (given in Eq. (1441)) is positive for 
L < 4 (s < s(g)) and becomes negative for L > 4 (s > s(q)). Hence, for 
all q > 1 the effective potential V(x) = /j,q + fi\x + [i2X q becomes unbounded 
from below when L exceeds 4. The solution of regime II is thus metastablc 
in the range s(q) < s < So(q) (4 < L < L ). Indeed, exactly as for q = 2, 
there exists a second solution for s > s(q) that becomes energetically more 
favorable (lower energy) for s > s 2 (<7). This is the onset of regime III. It 
occurs at s = S2 = s for very large N, more precisely at s — S2{q, N) ~ 

r -i2g/(2g-l) 

s(q) + yfqp(q ~ 1) s(q) /jv(8-i)/(2«-i) f or i arge but finite N, as we 

shall see. 

As the distribution of E g is highly peaked for large N, its mean value is 
given by the most probable value: (E g ) = N 1 ~ q s(q) where s(q) minimizes $(s). 
This minimum s = s(q) — ^^7^^ 4 9 (or equivalently L(s) = 4) is reached 

within regime II and $n(s(q)) = 0. For s close to s(q), $u(s) « ^"^f^ 

where a 2 is given in Eq. (|21l) . We conclude that the distribution of E g has a 
Gaussian behaviour around its average, as shown in Eq. (|20[) , from which we 
can read the variance (see Eq. (1211) ). For example, for q = 2, we have ct| = 4 
and VarE 2 

4.3 Regime III 

As s exceeds s(q), H2 becomes negative and the effective potential V(x) — 
Ho + \i\x + /J>2X q is not anymore bounded from below. The solution of regime II 
becomes metastable. The minimum of the potential at the origin still exists, as 
V(x) increases for small x, but it is a local minimum: V(x) reaches a maximum 
at x = x* > and then decreases to — oo (see Fig.[4jc)). Actually, for s > s(q), 
there exists another solution where one charge splits off the sea of the other 
(N — 1) charges that remain confined close to the origin (in the local minimum 
of V). The maximal eigenvalue (charge) becomes much larger than the other 
(see Fig. 12(c)). At some point s = S2(q,N) very close to s(q) for large N, this 
second solution becomes energetically more favorable than the solution of regime 
II : this is the onset of regime III. This phase transition occurs at s = S2(q,N) 
given in Eq. (|48j) . It is reminiscent of the real-space condensation phenomenon 
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observed in a class of lattice models for mass transport, where a single lattice 
site carries a thermodynamically large mass |38) . 



4.3.1 Regime III: summary of results 

We show in this section that there is an abrupt transition from regime II to III 
at s = s 2 (q, N) where: 

-,29/(25-1) 

/5/2( ff -l)«(<?) 

s 2 (q, N) « s(q) + J= ^-TTTToA &r large AT . (48) 



jV(g-l)/(2g-l) 

Here, s(q) the mean value of s given in Eq. (|19j) . The maximal eigenvalue t 
suddenly jumps from a value t ~ T/JV very close to the upper edge ( of the 
sea of eigenvalues to a value t « [s — ^(g)] 1 ^ 9 /N 1 ^^ much larger than the other 
eigenvalues (t 3> C) ( see Fig- 3 ( c ))- This is clearly shown by the good agreement 
between our predictions and numerical simulations in Fig. [7] for N — 500 and 
N = 1000. The consequence of this phase transition in the Coulomb gas is an 
abrupt change in the distribution of E g . More precisely, we show that for large 
N : 

P (Y, q = N 1 ' 11 s, N) « exp|-/37V 1+ 5 * 7J7 (s)| for s > s 2 (q, N) , (49) 
where 

*/7/(s) = . (50) 

The expression of the mean value s(q) is given in Eq. (jT9"]) . For example, for 
q = 2, this implies: 



p(E a = ^ J ^)«exp{- J 9iVi* J/ i(s)} with - . (51) 

The rate function $(N,s/N) defined by 

9 , , ( N 2 $ri(s) for s < s 2 , 

I iV + " *///(s) for s>s 2 , 

is continuous but its derivative is discontinuous at s — s 2 - for large N we have 
4^1 + ~ ^rl -/(2<?)- At the transition point s = s 2 , there is also a change of 

concavity of the curve: the rate function in regime II is convex ( d ?l 1 > for 
s < s 2 ) and has a minimum at s = s, whereas the rate function in regime III is 
concave ( d J g l" < for s > s 2 ). 

Figure [S] shows the transition from regime II to regime III for q = 2 and 
JV = 1000: analytical prediction for large N in Eq. (fSTj) compare well with 
Monte Carlo numerical simulations. 
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4.3.2 New saddle point 

We want to describe the regime where a single charge (the maximal eigenvalue) 
detaches from the continuum of the other charges. The assumption that all the 
eigenvalues are close to each other and can be described by a continuous density 
of states does not hold anymore. The saddle point must be slightly revised. 

We write X rnax = t and label the remaining (N — 1) eigenvalues by a continu- 
ous density p(A) = J2i^max ~~ ^»)- Physically, as the effective potential 
has a local minimum at the origin x — 0, we expect the optimal charge density 
p c to have a finite support over [0, (] with ( < t and p c (() — 0: while one charge 
(the maximal eigenvalue t) splits off the sea, the other charges (the sea) remain 
confined close to the origin (in the local minimum of V, see Fig. H](c)). 

In this regime, we do not rescale the density (and the energy) by assuming 
that A ~ 1/N. We want indeed to compute the pdf of E 9 = S for all S(q) < 
S < 1, where S(q) = N 1 ^ q s(q). The effective energy is now a function of both 
t and p: 

Es[p,t} = - (iV ; 1)2 f C [ C d\d\' p(X)p(X) ln|A-A'| 
1 Jo Jo 

- (N — 1) dX p(X) In \t - A| + p ^ jf C dX p(A) - l^J 
+ pi UN- l )f Q dX X P( X ) +t- lj 

+ p 2 (^{N -1) £ dXXi p(X)+ti - S^J . (53) 

The dominating configuration is described by the optimal charge density p c 
and the optimal value t c of t = X ma x such that: 



SE S 



5p 



= and 



p—p c ,t—t c dt 



= 0. (54) 

P=Pc ,t=t c 



Taking into account the normalization, we have indeed for large N: P (E q = S, N) 

! fvpidte-f>B\p't\ ~ ex P {-0 i E S [Pc t c ] - E [p*, t*])}, where E s [p, t] is given in 
Eq. (|53p and E [p, t] has the same expression as Es [p, t] but without the last 
term (the constraint ^ A| = S). The pair (p* ,t*) (resp. (p c ,t c )) minimizes 
E [p,t] (resp. Es [p, t]). In fact, the normalization is given by the saddle point 

energy evaluated at S — S (the mean value of S): E [p*,t*] = Es [p c ,t c 

(with S = 2/N for q = 2). We shall see that for large N, we have: 



E[p\t*}=Es[p c ,t c }\ s=5 ^N 2 (^f + ^j . (55) 

Formally, by analogy with regimes I and II, we can write: 

P (Eg = S, N) « exp {-f3N 2 <5>{N, S)} , (56) 
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where we define the rate function $ as 

$(AT, S) = (E S [p c , t c ] -E\p\ t* ]) /N 2 . (57) 

However, we shall see that the scaling of $ with N is different in regime III 
with respect to the regimes I and II. In regimes I and II, $ was independent 
of TV for large N: &(N,s/N) -> $/(s) (resp. <J>//(s)). In regime III, we shall 
see that: $(N,s/N) w ^i/j^/AT 1- * for large AT. 

For simplicity, we write £ instead of t c in the following. 

4.3.3 Case q = 2 

Following the same steps as for regime II, we find that the optimal charge 
density is explicitly given for q = 2 by: 

Pc{x) = ^hr)\[^ 

with A = 

C = J , where £ and t = t c satisfy: 



(a) 165 + N( 2 - 12C - ( 16t2 - 20t C + 5 C 2 ) = . ( 59 ) 

(b) (8t 2 - 8*C + C 2 ) 2 - 8(i - C) v 7 ^ - (8t - 2C - 2Nt( + N( 2 ) . (60) 

These equations can be solved numerically for every £2 = S. We can also find 
the solutions analytically for very large N. 

For S = with 2 < s < 9/4, there exist two solutions for the pair ((,t). 
The first solution is of the form t w £ with £ w 0(1/ N). This is exactly 
(to leading order in TV) the solution of regime II (see below, "first solution" ) . 
There is also a second solution, where t ^> (: the maximal eigenvalue becomes 
much larger than the other eigenvalues. More precisely, £ w 0(1/N) whereas 
t w 0(1/ VN) for 5 « 0(1 /N) (sec below, "second solution"). We shall see that 

the first solution (regime II) is valid up to a value s — s 2 ~ 2 + ^73- for large 
AT, whereas the solution with t ^> ( starts to dominate for s > s 2 (its energy 
becomes lower): this is regime III. 

For S > (s > |), there remains only one solution (the second one), where 
C = L/N and t > (. 

Note that in both cases, for large N (and for < S < 1), the upper bound 
( remains of the order ~ 0(1/N). We shall thus write C — ^ with L <~ 0(1). 
On the other hand, the maximal eigenvalue t scales from 0(1/ N) (as S —¥ 2/N) 
to 0(1) (as S -> 1"). 

Finally, we compute the saddle point energy as a function of £ = L/iV and i. 
As finite-size effects (large but finite N) are important in this regime, we keep 
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A + BX + 



C 
t-X 



(58) 



N(- 2 + 2y/t(t^O 



, B 



k (3C-4t) 



and 



all terms up to order O(N) in the saddle point energy, which gives: 



E s [p c ,t] = E((,t) 



(N-l) 2 



In 





- 2N In 




_4_ 




2 



9iV 2 



t-C 



6(1 + t 2 ) _ 5(N + t) 

lm m 11 



t 



8(t-C) 
5Nt~ 



C 



dn[t(t-C)] 



(61) 



where £ — ((s) and t = t c = t(s) are given by Eq. (59"]) and (fHU|) . 

The rate function is thus given by $(iV, S) = (E s [p c , t]-E[p*,t*\) /N 2 = 
(E [C, i] - E [p*,t*]) /N 2 with E [C, t] given in Eq. (|5T]>. 

Scaling S = s/N with s - 0(1) : first solution i f with C ~ 0(1/N) 
(regime II) 



For S 



N 



with s ~ 0(1) for large N, the solution of regime II still exists 



as long as s < 9/4 (where 9/4 = sq(2)). We recover this solution from the Eqs. 
(55} and (50]) with the scaling t = ^ and C = j} with T w L ~ O(l), i.e. the 
maximal eigenvalue i remains very close to the other eigenvalues (f ~ C f° r large 
AT). 

In this limit, equations ([59)) and (|60j) indeed give: 



(a) 16s + L 2 - 12L 
(6) {T-Lf /2 



0. 



L 5/2 



1 



8(6 -L) AT 



(62) 
(63) 



Equation (a) is the same as Eq. (14"5"]) of regime II. To leading order in N (order 
N 2 ), Eq. (HI]) reduces to: 



E s [ Pc ,t] = E(L,t) 



N 2 



In 



N 



5 ~ +N — + 8 



(64) 



Therefore, using Eq. (55]), we get $(AT, s/AQ = (£s[p c , i] - #s[Po * 
$(s) with $(s) = -|hi(f) 



6 



/N 2 = 

+ | = ( f>//(s). We recover the expression 

in Eq. (|4*2l of regime II. 

However, for S = s/N > 2/N there exists a second solution that becomes 
energetically more favorable at some point s 2 ~ 2 + . Therefore regime II 
is only valid for 5/4 < s < s 2 - 

Scaling S = s/A" with s ~ O(l) : second solution t C (regime III) 

For 5 = s/AT with s > 2, there exists a second solution where one eigenvalue 
(A max = t) becomes much larger than the others : t ^> (. In this limit, Eq. ([59]) 
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lnP(S 2 =f) 

f3N 2 

0.008 



0.006 - 
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2.0 2.1 12.2 2.3 2.4 2.5 



s(2) s 2 (2) S 

Figure 6: Distribution of E 2 : rate function $ = - In P [S 2 = s/N] /((3N 2 ) 
plotted against s for AT = 1000. Analytical results (solid line) are compared with 
data (red points) of numerical simulations (Monte Carlo, method 2, see section 
[5]) . Analytical results here are the rate functions expected in the limit of very 
large N: &u(s) in regime II (green solid line, see Eq. (H2"10 and $>(N, s/N) s» 
^in(s)/y/N in regime III (blue solid line, see Eq. (|51|) ). The transition between 
regimes II and III is abrupt, we can see the discontinuity of the derivative of 

the rate function. It occurs at S2(q — 2,JV) ~ 2 + jjtjj — 2 3^2/3^ ~ ^.18 f° r 
AT = 1000. 
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and (15(H) give for large AT: 



2 a r 4 
- and 



1 



3-s 



(65) 



= and 



For S — > 1, which implies s -> oo as JV -> oo, we find t R 
£ A (l — y/jf) s» ^ (1 — i) as also recovered in Eq. (fTUj) . 

We can expand the saddle point energy in Eq. (|6ip replacing i and £ by the 
expressions given in Eq. (|65|) for large AT. We obtain: 



E S \pc,t] 



v^2 



^3/2 + ^2 



In AT 1 



- ) - — In AT + O(iV) for large TV . 

(66) 

^3/2 _ 



Finally, we get N 2 <t>(N : s/N) = £ s/JV [p c ,*] - N 2 + ±) 
f In AT + 0(N) for large N (see Eq. l|gg]>) and the pdf of £ 2 is thus given for 
large N by: 

(67) 



P 



-/3AT 3 / 2 *ji X ( s ) 



where A^ 3/2 */7/(s) = AT 2 $(AT, s/JV), that is 



*///(s) = 



V7^2 IniV 



2VN 



O 



for large N . 



(68) 



The rate function has a very different behaviour for large N in regime II and 

py e~^ N *( s ), whereas in regime 

For large but finite iV and for 



III. In regime I and II, we have P (£2 = jj , N 
III we have P(E 2 = %,N) « e -/3JV 3 / 2 *„x( s ) 
s > 2 but very close to s = 2, we have A^ 3 / 2 *7 7/ (s) > N 2 $ n (s). Therefore 
the solution of regime II dominates close to s = 2. However, the solution of 
regime III becomes energetically more favorable at some point s 2 defined by 
N 3 / 2 ^ ln (s 2 ) = N 2 ^ II (s 2 ), that is 



2 4/3 2 5 / 3 lniV 



S2 



N l/3 3iV 2/3 



for large N . 



(69) 



At s = S2, there is an abrupt transition from regime II to III. The maximal 
eigenvalue t jumps from a value t s=s £ with T ~ 0(1) and i very close to £ 



to a value t 



•/N 



much larger than the other eigenvalues (t > £). The rate 



function is continuous but its derivative is discontinuous: A^ 2 rf ^ JJ 

(IS 



2273 ' 



whilst AT3/2rf£lil 

(IS 



iV 



5/3 



422/3 for large AT. At the transition point s = s%, 

there is also a change of concavity of the curve: the rate function in regime II 

s = 2, whereas 



is convex { a ^.l 1 > for all s < 9/4) and has a minimum at s 



the rate function in regime III is concave ( ^j" < for all s > 2). 
Scaling £ 2 = S ~ 0(1) and limit S — > 1 (unentangled state) 



2N 



In the far-right tail of the distribution £ 2 = S « 0(1) (S > s/JV, S < 1) 
and the maximal eigenvalue £ rj 0(1) whereas £ (and all the other eigenvalues) 
remain of order 0(1/N). In this limit, equations (|59"|) and (|60p become: 

5wt 2 and L«4(l-i) as i w O(l) . (70) 

The saddle point energy in Eq. (I5T1) reduces to: Eg [p c ,t] ~ — ^-ln(l — i) + 

iV 2 (^ + i)-AnnAT+O(A0asS«O(l)witht = VS. Using Eq. (|55|). we get 
an explicit expression for the rate function 3>(iV, S) = (i?s [p c , t c ) - E [p* , t*}) /N 2 
for large N: 

t(W , S) „ {gd*dz*ffl£±iB „ . 1 ,„ (l _ VS) s .,„<*, . (71) 
We conclude that 

P (E a =S,N)ki e-^ 2 * 111 ^ « (l - 2 for large iV, fixed 5 . (72) 

The difference of scaling with respect to regimes I and II comes from the scaling 
of S 2 : in regimes I (resp. II), we had §(N,s/N) — > $/(s) (resp. $//(s)) for 
large N, whereas here we have: $(N, S) — > $iu(S) for large N and fixed 
S « 0(1). As S = s/N with fixed s and large N, which corresponds to the 
limit S — > in this scaling, we find N 2 <P IH (S) « N 3 / 2 Vs/2 which is also the 
limit s oo of 7V 3/2 *///(s). The right tail (where S « 0(1/N)) and the 
far-right tail (where S « O(l)) of the distribution match smoothly. 

As £2 = S tends to its maximal value 1, the maximal eigenvalue t — > 1 and 
L — > 0. At S = 1, only one eigenvalue, the maximal one X m ax = is nonzero 
(and equal to one). This corresponds to an unentangled state (situation (i)). 
The probability of an unentangled state (i.e. £2 ~> 1) is thus vanishingly small 
for large N. 



4.3.4 General q > 1 

Using again Tricomi's theorem and imposing the constraints J p c = 1 and /? C (C)> 
we find that the optimal charge density for the N — 1 smallest eigenvalues is 
given by: 



Pc(A) 



?r(iV-l)V A 



C-A 



A + 5 2-F1 1,1 



3 A 



C 



t-X 



(73) 



where A = fix, B = p-ilqQ 1 1 T ^Y(q) anc ^ ^ = V t^c an< ^ * s a hypergeomet- 



ric function 2^1 (a, b, c, z) = Y^=o ^(c)^" S> with ( a )« = a(a + l)...(a + n - 1) 
denoting the raising factorial (Pochhammer symbol). The Lagrange multipliers 
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Nt = N\ r 
for E 




2.0 

s 



2.1 



12.2 

(N = 1000) 

S2 

(N = 500) 



2.3 



2.4 



2.5 



Figure 7: Maximal eigenvalue A max = t corresponding to a fixed value of the 
purity £2 = s/N plotted against s for different values of N. Analytical pre- 
dictions (solid lines) are compared with numerical simulations (points : Monte 
Carlo data, method 2 with density). The theory predicts for large N a sudden 
jump of t from a value t w £ = L(s)/N with L(s) = 2(3 — V9 — 4s) (within 
regime II, s < S2) to a much larger value i s» ^J^ 2 (regime III, s > s 2 )- We 
clearly see this jump in numerical simulations for N = 500 at S2 ~ 2.23 and 
N = 1000 at S2 ~ 2.18. For AT = 50, finite-size corrections to the large N 
asymptotics are considerable enough to smear the jump in t. Because of the 
choice of scaling on the plot, tN as a function of s, the plots of the maximal 
eigenvalue in regime II are expected to be the same for different ./V (for large 
N), whereas the plots for regime III differ by a factor \/N. 
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/ii and \i% are given by: 



Mi 



M2 



<77V C _2(g+l) + 



(9 - 1)C 2 

V^r( g + 2) 

C' +1 r(g+ l/2)g(g-l) 



t-C 



{(2q + 2)t-(2q+l)(} 



■i-N( + 



t-C 



{3C-44} 



where £ and t are solutions of the following system of equations: 



V ' v^r(g + l) \ 2 ^9+1 



-2i]} 



c 



2t(i - C) 



M2- 



C g r( g + i) 
t^Fr(g) 



(74) 



/3£+l 

c 



(76) 



with /^i = Mi(Cj*) an d M2 = M2(C>*) given in Eq. (|74|) . 

These equations can be solved analytically for large N and the solutions are 
qualitatively the same as for q = 2. 

For S = N 1 ^ s with s(q) < s < s (q) (where s (q) = 2 ^|+ 1 / + % (^)", 

see regime II), there exist two different solutions for the pair ((,t). The first 
solution is of the form t pa £ with £ ps 0(1 /N). This is exactly (to leading 
order in AT) the solution of regime II (see below, "first solution"). There is 
also a second solution with t £> more precisely £ = L/N with L ~ 0(1) 
and i » CKl/iV 1 - 1 /-?) for S ps A^- 9 s, and C « 0(1/AT) (see below, "second 
solution"). For s close to s(q), the first solution dominates (regime II), but at 
some point s — s 2 (q,N) > s(q) given in Eq. (|80|) . the second solution, with 
t ^> Cj starts to dominate (its energy becomes lower): this is regime III. 

For S > N 1 ~ q s (q), i.e. s > So, only the second solution remains: the 
upper bound of the density support scales as ( — L/N with L ~ O(l) while the 
maximal eigenvalue is much larger than all other eigenvalues: t 3> C- 

In both cases (as for q = 2), for large N the upper bound ( remains of order 
- 0(1 /N) (C ~ A typ ). We shall thus write C = j? with L ~ 0(1). On the 
other hand (as for q = 2), the maximal eigenvalue t scales from 0(1/ N) (as 
S -> JV x -«5(g)) to 0(1) (as 5 -> 1"). 



Scaling S = A^ 1 - 
0(1/AT) (regime II) 



9 s with s ~ 0(1) : first solution t ps £ with £ 



For 5* = A^ 1-9 s with s ~ O(l) for large AT, the solution of regime II still 
exists as long as s < sq (q) . We recover this solution from the Eq. (|75j) and ([TBI 
with the scaling t = jt and ( = -k with T pa L ~ O(l), where the maximal 
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eigenvalue t remains very close to the other eigenvalues (t « £ for large TV), 
it does not play a special role. Using Eq. (|55|) . we finally get <&(iV, s/JV) = 

(Es[p c ,t] — Es[p c ,t] ) /N 2 = 3>jj(s), the same expression as in Eq. (|37|) of 

regime II. 

However, for s > s(q) there exists a second solution that becomes energeti- 
cally more favorable at some point s%(q,N). Therefore regime II is only valid 
for si < s < S2- 

Scaling S = iV 1 ^ 9 s with s ~ 0(1) : second solution t > ( (regime III) 

For S = N 1 ~ q s with s > s(q), there exists a second solution where one 
eigenvalue (A max = t) becomes much larger than the other eigenvalues : t S> C- 
In this limit, Eq. ([75} and (|75|) give for large AT: 



fa)] 



1/9 



and 



s-«(g)(l + <?)/2 



1 



%)] 



1-1/9 [ ATi-i/9 



(77) 
51/9 



For 5—^1, which implies s -> 00 as iV -> 00, we find £ « s 1 / 9 A^ 1 / 9 
and C «£(!-*)■ 

We can compute the saddle point energy in this limit replacing t and £ by 
the expressions given in Eq. (f77f for large N. Finally, we get N 2 Q(N,s/N) = 

E s [ Pc , t] - N 2 + |) « iV 1+ « [s ~ g( 2 g)ll/g for large AT (see Eq. |55j) and the 
pdf of Eg is thus given for large N by: 



where 



P (Eg = N 1 -* s, N) a exp {-£AT 1+ 5 #/r,(s)} , 
[ S -%)] 1/9 



*///(s) 



for large N . 



(78) 



(79) 



The solution of regime III becomes energetically more favorable, that is 



N 1+1 "^ ln (s) < iV 2 $//(s), at some point s 2 (q,N) defined by N 1+ ? *///(s 2 ) 
N 2 $ n (s 2 ). Therefore 



i+i 



s 2 (q,N)^s(q) + 



^/2 (ff -!)«(«) 



29/(29-1) 



AT(9-l)/(29-l) 



for large N . 



(80) 



At s = S2, there is an abrupt transition from regime II to III. The maximal 
eigenvalue t jumps from a value t m with T ~ O(l) and i very close to £ to 

a value t « ^jvi^i//" mucn larger than the other eigenvalues (t> C) ■ The rate 
function Q(N,s/N) given by 



N 2 $(AT, s/AT) 



A' 2 <frj/(s) for s < s 2 , 
N 1+ * *///(s) for s > s 2 , 



(81) 
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is continuous but its derivative is discontinuous. For large N, we have indeed 

' -2q 

T 2 d<& I 



^ ds 



« N m {( q - 1)^/2 s( q )Y q -\ whilst f | s + a f L-/(2«). At 
the transition point s = S2, there is also a change of concavity of the curve: the 
rate function in regime II is convex ( d ^ I 1 > 0) and has a minimum at s = s, 
whereas the rate function in regime III is concave ( d < 0). 



5 Distribution of the Renyi entropy S q 

In section [H we have computed the full distribution of E 9 = J2i=i tf, f° r large 
N. A simple change of variable gives the distribution of the Renyi entropy 



S a 



1-9 



In [£ q ]. The scaling E 9 = iV 1 q s for large N implies 



S q = In N — ^Mj- . This means that typical values of S q will be of order S q s» 
lnAT — z with z s» O(l) for large N. The parameter z = ^j- is nonnegative 
and its minimum z — corresponds to = IniV, which corresponds to the 
maximally entangled state. 

The distribution of the entropy is thus given for large N by: 



P(S q = hxN-z) 



exp{-/3iV 2 <f>x(g)} 



for < z < zi(q) 



exp {~/3N 2 (j)ii(z)} for zi(q) < z < z 2 (q) , 



exp 



{-/3iV 1+ i fc,(z)} for z>z 2 (g). 



The three regimes are the same as for E„. The rate functions </>/, 



(82) 
6// and 



ipiu are simply obtained from the rate functions 3>jj and ^j// for the 
distribution of E g (see Eq. ([18])') by the change of variable s = exp [(g — l)z], 
e.g. 0/(z) = $/ (e^- 1 ) 2 ). Explicit expressions of the functions $/ and $// are 
given in Eq. (13"5]) and (H2J) for g = 2, and in Eq. (|T7]) for a general q > 1; 
an explicit expression of ^ m is given in Eq. (|50[) for a general g > 1 (and in 
Eq. AID for g = 2). 

The critical points are given by 



lnsi (?) lns 2 q,7V 

ziW = T and z n1> N ) = 1 

q-l q-1 



(83) 



where s± and s 2 are the critical points for E 9 (see Eqs. (|22|) and ([23])). 

The distribution of the entropy S q has the same qualitative behaviour as 
that of E g : it is a highly peaked distribution with Gaussian behaviour around 
the mean value and non-Gaussian tails. Again, the average value of S q coincides 
with the most probable value for large N, (S q ) s» ln N — z(q) where z(q) is the 
minimum of 4>n- 



(S q ) fa In N - z(q) with z(q) = 



1 



ln s(q) 
q-l ~ 9-1 



In 



r(g + l/2) A" 
T(q + 2) ^ 



■ (84) 
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The rate function 4>n(z) has a quadratic behaviour around z = z(q): cj>u(z) ~ 
- — q " . Therefore, the distribution of the entropy S q has a Gaussian behaviour 
around its average: 

P(S q =hxN-z) saexp|-/37V 2 (z - Z ^ \ for z « 2(g) , (85) 

which gives the variance of the distribution: 



Var5 9~ forlar S e7V - ( 86 ) 



5.1 Limit g — >■ 1 + : von Neumann entropy 

As q — > 1 + , the Renyi entropy S q tends to the von Neumann entropy Svn = 
— J^i \ hiAj. The limit q 1 is singular for the distribution of Y> q : because 
of the constraint Si = Aj = 1, the distribution tends to a Dirac-(5 function. 
The variance tends to zero (a 2 — > 0) and the mean value s(q) as well as the 
critical point si(q) and S2(q) tend to 1. However, due to the factor 1/(1 — q) 
in the definition of S q , the limit q — > 1 is not at all singular for the entropy S q . 
Taking this limit only requires to be careful. For Svn (as for S q for q > 1), 
there are three regimes in the distribution: 



P(S VN =UN-z) 



exp {-/3N 2 4>i{z)} for < z < z x , 
exp {-/3N 2 4>ii{z)} for zi < z < z 2 , (87) 
exp 0///(^)| for z > z 2 , 



where 0// and 0/// are respectively given in Eqs. (|9T))) and (|9~4")l . For g — > 1, we 
get: 2(g) — — >• 1/2 (where 2(g) is given in Eq. ([Ml) ). We thus recover 

the already known mean value of the von Neumann entropy (see [3]) in the case 
c = 1 (M re AT): 

(Svn) « IniV - - for large N . (88) 

The critical points separating the three regimes are given by (limit q — > 1 in 
Eqs. ([13} and 

2 3 1 
Z! = - - In - w 0.26 and z 2 ~ 2 = - . (89) 

We easily obtain the expression of the rate function <pn in regime II by 
taking the limit q — > 1. We get: 

^(*) = 4 ln (z) + A _ z + 1 ' (90) 
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where L = L(z) is the solution of (limit q — > 1 in Eq. (|4"B"|) ) 



For large AT, the mean value corresponds to the minimum of tfijj. The quadratic 
approximation of (fin around this minimum zkz gives the Gaussian behaviour 
of the pdf of 5*vn around its average and thus the variance in the large N limit: 

(Svn) ~\nN-z with z = - and VarS V N « ■ (92) 

The limit q — > 1 for the regime III is a bit more subtle. We would expect 
the rate function to be of the form N 2 ip HI (z), but ifim = \Pjxr(e( ?_1 )*) (in Eq. 
([5D|) ) vanishes as g — > 1. The rate function actually scales as N 2 /\nN (rather 
than N 2 as one could naively expect). This can be shown by a more detailed 
analysis of the equations (|75|) and (|76|) in the limit q — > 1 . The solution £ 3> C 
is actually given for q — > 1 by: 

z - 1/2 4 / 1 - z\ 

In AT s TV V In AT / v ; 

The saddle point energy can be computed in this limit. We finally find: 

iV 2 1 
- lnP(SvN = In AT -z) 1/2); <fiiu{z)=z-~. (94) 

5.2 Limit q — > oo : maximal eigenvalue 

As g — > oo the Renyi entropy S q tends to — In A max where A max is the maximal 

eigenvalue. Again, the limit is singular for the distribution of S g but not for S q . 

There are the same three regimes in the distribution of A max for large A" as in 

the distribution of the Renyi entropy. 

For large N, the typical scaling is S q « In A~ — z, thus — In A max w InJV-z or 

A max ~ jj- Setting t = e 2 , we have A max = t/N. In particular, the mean value 

_ _ _ _ i 

is given by t/N where t = lim^oo exp(z(q)) = lim q -nx> [s(q)\ "- 1 = 4, implying 

(A max ) « 1 . (95) 
i 

The first critical point is t\ — lim^oo [si(q)] q ^ 1 = 4/3. The second critical 
point is £2 = t = 4. The three regimes in the distribution of the maximal 
eigenvalue are the following: 



t 

N 



e -PN 2 X r{t) for i <t < 4/3 (reg. I) , 
e -0N\ ri (t) for 4/3 < i < 4 ( reg . n) , (96) 
e -PN XlII (t) for t > 4 ( reg _ ni )_ 
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The rate functions can be explicitely computed. The rate function in regime I 
is given by: 

Xi(t) = ~Ht-l) forl<t<4/3. (97) 
In regime II, we hnd: 

Xu{t)=^^^-\\n(^j+\ for4/3<i<4. (98) 

Finally, in regime III the maximal eigenvalue detaches from the sea of the other 
eigenvalues and we get: 



Xm(t) = — ^ 2 ^ -2rn(V7+VT^~4) + 21n2 for t > 4 . (99) 

Again, at the first critical point t\ = 4/3, the rate function x is continuous and 
twice differentiable, but its third derivative is discontinuous: = —27 but 

d X" = —999/64. The average value t = 4 is the minimum of Xn- At the second 
critical point i 2 = 4, the rate function is continuous but not differentiable. 

Exactly as we did for T, q , we can also consider the regime where A max = T 
(T » i/AT): the far-right tail of the distribution. We find: 

P (A max = T) « e-? N2 x+m X+(T ) = -I ln(l - T) for < T < 1 , (100) 

which matches smoothly regime III. We have indeed: Nxin{t) ~ AT| as t — > oo 
and A f2 x+(t) ~ N 2 ^ m asT ^0 with T = t/iV. 

Ideas of proof 

Regimes II and III can be derived by taking carefully the limit q — >• oo 
(directly in the expression of the rate function for regime II but more carefully 
for regime III). The distribution of A max can also be computed directly (without 
taking the limit q — >• oo). This gives the same results for regimes II and III 
and gives also an explicit expression for regime I (where the rate function is 
not explicitely known for a general q > 1). We can actually calculate the 
cumulative distribution Prob (A max < Z) by the same Coulomb gas method as 
before. This is indeed easier to compute because the probability that A max < Z 
is the probability that all the eigenvalues \ are smaller than Z. We can thus 
compute this probability with the Coulomb gas method, with a continuous 
density p(x) = 1/N S(x — XiN) and with the constraint that no eigenvalue 
exceeds Z: 

P(A max < Z) oc / Vpe-P"** 1 *™ . (101) 
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The energy reads 



Ez[p] 



— / p(x)p(x') In \x — x'\ dx dx 1 + /i / p(x)d. 
2 Jo Jo \Jo 



1 [ z [ z 




■X 



1 



) 



i / xp(x)dx — 1 




(102) 



where the Lagrange multipliers po and p\ enforce the two constraints J p = 
1 (normalization of the density) and J xp = 1 (unit sum of the eigenvalues: 
J2i = !)■ The saddle point method gives: 



P(A : 



< Z) on e 



-l3N 2 E z [ Pc \ 



(103) 



i max 



where p c minimizes the effective energy Ez- This yields regimes I and II. 
Exactly as for S q , in regime III, the maximal eigenvalue detaches from the sea 
of the other charges (eigenvalues), it must be taken into account separately from 
the continuous density of the other eigenvalues. 

In regime I, the optimal charge density has a finite support [Li,^] and 
vanishes at L x 2 (exactly as for S g ). We get the rate function \i m Eq. (jHTJ). 

In regime II, the optimal charge density has a finite support ]0, L], vanishes 
at L but diverges at the origin with a square root divergence (exactly as for 
S 9 ). We get the rate function xii m Eq. (f^5|) . This expression can also be 
obtained by taking the limit q — > oo of the expression in Eq. (j47]) of valid 
for a general q (for E g ). 

In regime III, the maximal eigenvalue is much larger than the others and we 
get xiii m Eq. (|99p . The limit q — > oo in the rate function ipm for a general q 
gives: ipm — > t/2. This is actually equal to Xin(t) only in the limit t — > oo, 
but not for all t > 4. For q > 1, regime III is characterized by t ~ T/N 1 ^^ ^> ( 
as C ~ L/N, which becomes t « T/N > £ in the limit q — > oo. The maximal 
eigenvalue is larger than the other eigenvalues, but not much larger. We cannot 
anymore assume t 3> C m t ne computation of the energy. We must compute 
carefully the energy Eg [p c ,t] in this limit. For this computation, we use the 
complete expression of Eg: for q = 2, this expression was given in Eq. (1611) : for 
a general q, we have a similar but more complicated expression. We use this 
expression in the limit where t and £ are both of order one (with t > £) and 
where q — > oo. We finally get xm(t) as given in Eq. 



5.2.1 Typical fluctuations around the average: Tracy- Widom distri- 



We have seen that the average value of the maximal eigenvalue, in the large N 
limit, is given by (A max ) s» 4/JV. Of course, A max fluctuates around this average 
from sample to sample. The Coulomb gas method presented in this subsection 
captures fluctuations ~ 0(1/N) around this mean, i.e., large fluctuations that 



bution 
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are of the same order of magnitude as the mean itself. We have seen that the 
probability of such large ~ 0(1 /N) fluctuations is very small, indicating that 
they are rare atypical fluctuations. The typical fluctuations around the mean 
occur at a much finer scale around this mean which is not captured by the 
Coulomb gas method. 

To compute the distribution of such typical fluctuations, we start from the 
joint distribution in ([5]). The cumulative probability of the maximum can be 
written as the multiple integral 



pZ 

P(X mscx <Z)(X / .../ P(\ 1 ,\2,...,\N)d\idX 2 --.d\ N 
Jo Jo 



(104) 



Next we can replace the delta function 6 Q^iLi A« — l) by its integral represen- 
tation: 5(x) = (l/2ni) J dpe px where the integral runs over the imaginary axis. 
This gives, for M = N, 



P (A max < Z) OC 



dp 



2ni 



N 



N 



H</\ « ^ A II A ' 1[A, A, . 

i°> z ] Li=l J i—1 i<j 

(105) 

Rescaling \ — > (f3/2p)Xi, one can recast the integral as 



P (A max < Z) 



d P C P v ~PN 2 /2 

2m 1 



[0,2pZ/P] 



N 



N 



3 \ 

e _ "2 £->i=i A 



i—i i<j 

(106) 

The integral over A^'s is just proportional to the cumulative distribution of 
the maximum of the Wishart matrix, i.e., the p^ shart (A max < 2pZ/0). This 
latter quantity, in the large N limit, is known [311130] to converge to a limiting 
distribution known as the Tracy- Widom distribution |41j . i.e, 



P 



Wishart 



(A max <V) Fp 



(y - 4JV) 

24/3JV1/3 



(107) 



where Fp(x) satisfies a nonlinear differential equation |41j . Using this result in 
P06p . we get, in the large N limit, 



P (A r „ax < Z) (X 



dP_ e P-§N 2 log(p) p 

2-ni p 



^Z-AN 

2 4 /3ATl/3 



(108) 



The integral over p can now be evaluated via the saddle point method. To 
leading order for large N, one can show that the saddle point occurs at p* = 
0N 2 /2 that just minimise the exponential factor e p ~? N log ( p ) . Hence, to leading 
order in large N , we obtain our main result 



P (Amax < Z) « Fp 



Z-A/N 

2 4/3jV-5/3 



(109) 
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This shows that A max in our problem typically fluctuates on a scale 0(N 5 / 3 ) 
around its average 4/ AT, 

typical A max = -i + 2 4 / 3 N- 5 / 3 X p , (110) 

where the distribution of the random variable xp is the Tracy- Widom probability 
density function gg(x) = dFp(x)/dx. Around the mean value we have then 

P (A max = £) * N^gp (2-^NV\t 4)) . (Ill) 



Matching between the tails of the Tracy- Widom distribution and the 
large deviation rate functions 

For Gaussian and Wishart matrices, it has been recently demonstrated [26l 
|27|, [28] that the Tracy- Widom density describing the probability of typical fluc- 
tuations of the largest eigenvalue matches smoothly, near its tails, with the left 
and right rate functions that describe the probabilty of atypical large fluctu- 
ations. It would be interesting to see if the same matching happens in our 
problem as well. Indeed, we find that the tails of the Tracy- Widom distribution 
match smoothly to our previously obtained rate functions. 

For the left tail of the Tracy- Widom density, it is known [JT] that gp(x) ~ 

cxp{-^|x| 3 } for x -> -co. Therefore P (A max = ~ exp { -/?AT 2 }. 

On the other hand, for the rate function to the left of the mean describing large 
fluctuations of ~ 0(1 /N) is given in (j98|) . Taking the limit t — > 4 - , we find 

Xn(t) ~ — ^384 thus matching smoothly with the left tail of the Tracy- Widom 
density. 

For the right tail, one knows |JT] gp(x) ~ exp | — ^a; 3 / 2 | for x — > +oo. 

Therefore P (A max = -£) ~ exp {-/3AT ttdj)— j. On the other hand, the rate 
function describing large fluctuations of order ~ 0(1/ N) to the right of the 
mean is given in ((99]) . Expanding to leading order for t — > 4 + , wc get: xni(t) ~ 

i — J — which clearly matches smoothly to the right tail of the Tracy- Widom 
density. 



6 Numerical simulations 

To verify the analytical predictions derived in the preceding sections, we simu- 
lated the joint distribution of eigenvalues in Eq. (0: 

P(X 1 ,...,X N ) = 5M,^(l>-l) f[A/ (M - Ar+1) - 1 niAi-A i |' 3 

\ i / i=l i<j 

= B M , N 5(^^-1) e~? E ^, (112) 
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where the effective energy E [Aj] is given by Eq. (fTTj) . We sampled this proba- 
bility distribution using a Monte Carlo Metropolis algorithm (see |42j). 



6.1 Standard Metropolis algorithm 



We start with an initial configuration of the A^'s satisfying = 1 and 

Aj > for all i. At each step, a small modification {A^} — > {A-} is proposed 
in the configuration space. In our algorithm, the proposed move consists of 
picking at random a pair (Aj, Afc) (with j ^ k) and proposing to modify them as 
(Xj, Afc) — > (Aj+e, Afe— e), which naturally conserves the sum of the eigenvalues, 
e is a real number drawn from a Gaussian distribution with mean zero and with 
a variance that is set to achieve an average rejection rate 1/2. 

The move is rejected if one of the eigenvalues becomes negative. Otherwise, 
the move is accepted with the standard probability 

P [^-^'n) A = min (^[{xryj-E^}]) A (U3) 
P(Ai, Xn) / \ > 

and rejected with probability 1 — -p. This dynamics enforces detailed balance 
and ensures that at long times the algorithm reaches thermal equilibrium (at 
inverse "temperature" (3) with the correct Boltzmann weight e~^ E ^*^ and 
with J2i Ai = 1. 

At long times (from about 10 6 steps in our case), the Metropolis algorithm 
thus generates samples of {Ai} drawn from the joint distribution in Eq. (|1 12[) . 
We can then start to compute some functions of the A,'s, e.g. the purity S 2 
V Af , and construct histograms, e.g. for the density, the purity, etc.. 

However, as the distribution of the purity (as well as the one of the eigen- 
values) is highly peaked around its average, a standard Metropolis algorithm 
does not allow to explore in a "reasonable" time a wide range of values of the 
purity. The probability to reach a value £2 = s/N decreases rapidly with N 
as e~^ N *M where $(s) is a positive constant (for s different from the mean 
value : s^s). Therefore, we modified the algorithm in order to explore the full 
distribution of the purity and to compare it with our analytical predictions. 



6.2 Method 1 : Conditional probabilities 

It is difficult to reach large values E 2 = s/N of the purity (s > s). The idea 
is thus to force the algorithm to explore the region s > s c for different values 
of s c . We thus add in the algorithm the constraint s > s c . More precisely, we 
start with an initial configuration that, in addition to A^ = 1 and A^ > for 
all i, satisfies also ^ A| > s c /N. At each step, the proposed move is rejected if 
J2i Xf < s c /N. If J2i Af > s c /N, then the move is accepted or rejected exactly 
with the same Metropolis rules as before. Because of the new constraint s > s c , 
the moves are rejected much more often than before. Therefore the variance of 
the Gaussian distribution P(e) has to be taken smaller to achieve a rejection 
rate 1/2. 
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We run the program for several values of s c (about 20 different values) and we 
construct a histogram of the purity for each value s c . This gives the conditional 
probability distribution P (£ 2 = ■§ , |^2 > jtf)- Again, as the distribution of the 
purity is highly peaked, the algorithm can only explore a very small range of 
values of s - even for a large running time (about 10 s steps). The difference 
with the previous algorithm is that we can now explore small regions of the 
form s c < s < s c + rj for every s c , whereas before we could only explore the 
neighbourhood of the mean value s. 

The distribution of the purity is given by 

P(S 2 = ^)=P(S 2 = ^|E 2 >|)*P(S 2 >|) (for Sc < S ). (114) 
Therefore the rate function reads: 



1 



(3N 2 



h,/ME2 = ^ S2 ^!) +lnP ( S2 ^ 



(115) 



The histogram constructed by the algorithm with the constraint s > s c is the 
rate function $ Sc (s) = —jj^ InP (S 2 = ^|S 2 > 3> Sc (s) differs from the 
exact rate function $(s) by an additive constant that depends on s c . In order 
to get rid of this constant, we construct from the histogram giving $ Sc (s) the 
derivative rf ^ c . This derivative is equal to d j and the constants disap- 
pear. We can now compare numerical data with the derivative of the analytical 
expression for the rate function $(s). 

We can also come back to 3>(s) from its derivative using an interpolation 
of the data for the derivative and a numerical integration of the interpolation. 
This allows to compare directly the numerical results with the theoretical rate 
function $(s). 

We can follow the same steps to explore the region on the left of the mean 
value s < s by adding in the simulations the condition Xf < 4$ (instead of 
Yli > ]y) f° r several values of s c < s. 

We typically run the simulations for N — 50 and 10 8 iterations. As figure 
[5] shows, numerical data and analytical predictions agree very well for regimes 
I and II (rate functions given in Eqs. (|38|) and (|42j)). For regime III, finite- 
size effects are important and agreement holds for large but finite N analytical 
formulae (taking as rate function the expression of the energy in Eq. (|6ip with 
t and £ numerical solutions of the system of equations ([5^]) and The 
agreement would degrade for the asymptotic rate function giving only the dom- 
inant term for very large N ( Eq. ([BT]) ). Finite-size effects are also important 
for the transition between regimes II and III. Large- TV data are crucial to sec 
clearly this abrupt transition with a sudden jump of the maximal eigenvalue. 
For N = 50, the transition appears indeed to be smoothed out. This obser- 
vation can be rationalized as follows. At the transition (s = s 2 ), the maximal 
eigenvalue t is expected to jump for large N from a value ~ to a much larger 
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yet for N = 50 we have > \ hjr for all s < 9/4. We thus 



value - 

conclude that no jump can be seen at N — 50 and much larger N are needed. 
Adapting the simulation method to cope with this challenge is the subject of 
the next subsection. 



6.3 Method 2 : Simulation of the density of eigenvalues 
(and conditional probabilities) 

We want to be able to run simulation for very large values of N. The idea 
is to simulate the density p(X) = ^\ — ^») rather than the eigenvalues 
themselves. In the previous scheme, a configuration was made of N variables, 
the N eigenvalues. In the new code, we have k + 2 <C N variables: 

(1) the maximal eigenvalue t. 

(2) the upper bound of the density support £ (( < t). 

(3) the value of the density at each point Xi = (for < i < k). 

We must enforce the condition p(Q = 0, i.e. p(xk) — by definition of the 
upper bound £ of the density support. The idea is to replace the real density by 
a linear approximation of the density defined by its value at Xi for < i < k. 

These k + 2 variables describing the maximal eigenvalue and the density 
of the other eigenvalues simulate configurations with N ^> k eigenvalues, for 
example N = 1000 with k = 50. The number of eigenvalues N appears in the 
expression of the energy (and in the constraints). With this new code, we can 
now simulate configurations with many eigenvalues in a reasonable time. 

The algorithm 

From the analytical calculations, we expect that the density diverges when 
A — > + as p(\) <~ In order to get a better approximation in our code, 

we choose to discretize a regularized form of the density p(X) = VAp(A). Our 
(k + 2) variables are thus: 

(1) the maximal eigenvalue t. 

(2) the upper bound of the regularized density support ( (( < t), which is 
the same as the upper bound of the density support. 

(3) the value of the regularized density at each point Xj = ^ (for < i < k): 
z t = p(x t ). 

In the Monte Carlo simulation, we compute the energy as well as the con- 
straints = 1; etc by using a linear interpolation of the regularized 
density p(X) : 

p(X) = z t + -^ii —(X — Xi) for A e [xi,x i+1 [, (116) 

Xj-f-i Xi 

with Zi = p(xi) (in particular Zk = 0). Integrals such as J dXXp(X) are com- 
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putcd using the linear interpolation as : 



Jo 



d\p(X)X « ± 2 z + |> {(i + 1)5 + (i - 1)1 - 2,1} 
L *=i 



• (H7) 



There are two constraints for the density: the normalization J p = 1 and 
the unit sum of the eigenvalues £ + (N — 1) / Xp = 1. We start from an initial 
configuration satisfying these constraints : for example, we can take for the ini- 
tial p a density of the form of the (normalized) average density p(X) = ^-\J ^j^- 

and fix t with the unit sum constraint t = —(N — 1) J Xp + 1. Initially, we also 
choose ( not too large such that the condition ^ Xf > s c /N is satisfied (for a 
fixed value of s c ), exactly as in the code with conditional probabilities. 

At each step, we propose a move in the configuration space (our k + 2 
variables) that naturally enforces the two constraints J p = 1 and t + (N — 
1) f Xp = 1 (unit sum). More precisely, at each step we choose randomly three 
integers between and k + 1 : i\ < i 2 < i 3 . 

• If is < k (case 1), we propose a move (z^, Zi 2 , Z{ 3 ) — > (z,i + aie,Z{2 + 
&2€,Zi3 + a^e), where e is drawn from a Gaussian distribution with zero 
mean and a variance adjusted to have the standard rejection rate 1/2 
at the end. a\, 0,2 and C13 are constants that are chosen such that the 
constraints f p = 1 and t + (N — 1) J Xp = 1 (unit sum of eigenvalues) are 
satisfied: 

\h + 1) 3/2 + (*3 - 1) 3/2 - 2*3 /2 l \(i2 + 1) 5/2 + (*2 - 1) 5/2 - 



a 1 



(i 2 + 1) 3/2 + (i 2 - 1) 3/2 - 2z 3/2 l [(i 3 + I) 5 / 2 + (i 3 - 1) 5/2 - 2* 



.5/2 



a 2 and 03 are obtained from ai by cyclic permutation of i\, i 2 and Z3. 

If ii < 12 < «3 = k (case 2) , we propose a move ((, , 2;^ ) — ► (^ + e, ^ji + 
ei, z%2 + £2) where e is drawn from a Gaussian distribution with zero mean 
and a variance adjusted to have the standard rejection rate 1/2 at the end 
(different from the variance of case 1), and where ei and £2 are functions 
of e, i\ and 12 fixed by the two constraints (/ p — 1 and unit sum). 

If ii < %2 < k and i 3 = fc + 1 (case 3), we propose a move (i, 2^ , 0j 2 ) — ^ 
(t + e, zn + e\,Zi2 + £2)7 where, exactly as in case 2, e is drawn from a 
Gaussian distribution, and t\ and €2 are functions of e, ii and ?2 fixed by 
the two constraints ( J p = 1 and unit sum). 

• If ii < Z2 = fc and 13 = fc + 1 (case 4), we propose a move (£, Zj-^t) — > 
(( + e,zn + ei,i + dt), where e is drawn from a Gaussian distribution 
(same as in case 2), and ei and dt are functions of e and i\ fixed by the 
two constraints (J p = 1 and unit sum) . 

Then, if ( > t, if C < 0, if z t < or if ^ A 4 2 < s c /iV, that is (N - 1) / A 2 p + 
t 2 < s c /N, the move is rejected. Otherwise we compute the energy of the new 
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configuration E ncw and accept the move with the usual Metropolis probability 
p = min (e~^ Enmr ~ B \ l) (and reject it with probability 1 — p). 

Direct inspection of the previous rules shows that detailed balance is sat- 
ified. Therefore, after a large number of iterations, thermal equilibrium with 
the appropriate Boltzmann weight is reached and we can start to construct his- 
tograms of the density and the purity. We verified that for N = 50 (simulated 
with k + 2 variables, where k = 20) we recover the results of the direct Monte 
Carlo (where we simulate directly the eigenvalues). For N = 500 and N = 1000 
(with k = 50), we get very interesting results that can be used to test the 
large- N analytical predictions (see Eqs. (|38p and (|4"2"j) for regimes I and II and 
Eq. (fSTj) for regime III) : figure [5] shows the good agreement between theory 
and numerical simulations with this second method, for the distribution of the 
purity S 2 = J2i with N = 1000. As figure [7] shows, we can really see the 
abrupt jump of the maximal eigenvalue and the change of behaviour of the rate 
function (discontinuous derivative) , which is expected at the transition between 
regime II and regime III for very large TV. 

The simulations also provide solid support to the fact that a single eigen- 
value detaches from the sea in regime III. One might indeed wonder whether 
configurations with multiple charges detaching from the sea could be more fa- 
vorable. This was ruled out by measuring the area of the rightmost "bump" in 
the density of charges (see Fig. [3]) and verifying that it corresponds to a single 
charge. This fact is also intuitively rationalized as follows. Let us consider con- 
figurations with two charges, Ai and A2 (Ai > A2), detaching from the sea. As in 
Eq. ([53]) . we require Af + A2 = t q and we consider the quantity C = 1 — Ai — A2, 
which quantifies the compression of the sea of charges and would replace 1 — t in 
the m constraint in Eq. ([53]). The smaller is C, the stronger is the compression 
of the sea (with the other constraints remaining the same). Since the charges 
repel each other, the energy of the configuration is expected to increase as C 
gets smaller. An elementary calculation shows that, due to the convexity of A 9 
for q > 1, C is minimum when Ai = A 2 = 2~ 1 / q t while its maximum (minimum 
energy) is attained at the boundary Ai = t, A2 = 0, corresponding indeed to a 
single charge detaching from the sea. 

7 Conclusion 

In this paper, by using a Coulomb gas method, we have computed the distribu- 
tion of the Renyi entropy S q for q > 1 for a random pure state in a large bipartite 
quantum system, i.e. with a large dimension N of the smaller subsystem. We 
have showed that there are three regimes in the distribution P (S q = In N — z) 
that are a direct consequence of two phase transitions in the associated Coulomb 
gas. 

(i) Regime I corresponds to the left tail of the distribution (0 < z < z\ (q)). 
In this phase, the effective potential seen by the Coulomb charges has a minimum 
at a nonzero point. The charge density has a finite support over [Li,L2[ (and 
vanishes at L\ and L2), the charges accumulate around the minimum of the 
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potential. 

(ii) Regime II describes the central part of the distribution (zi(q) < z < 
Z2(q)), and in particular the vicinity of the mean value z(q). At the transition 
between regimes I and II, the third derivative of the rate function (logarithm of 
the distribution) is discontinuous. In this phase, the charges concentrate around 
the origin, the charge density has a finite support over [0,L] with a square-root 
divergence at the origin. Close to the mean value of S q , the distribution is 
Gaussian. 

(iii) Regime III describes the right tail of the distribution (z > Za(<?))) 
corresponding to a more and more unentangled state. In this phase, one charge 
splits off the sea of the other charges. The transition between regimes II and 
III is abrupt with a sudden jump of the rightmost charge (largest eigenvalue). 
There is thus a discontinuity of the derivative of the rate function and the scaling 
with A changes at this point. 

A by-product of our results is the fact that, although the average entropy is 
close to its maximal value In N, the probability of a maximally entangled state 
is actually very small. The probability density function of the entropy indeed 
vanishes at z — (far left tail), i.e. at S q — In A, which is the maximally entan- 
gled situation. Similar properties and three different regimes are also obtained 
in the limit q — > 1, which gives us the distribution of the von Neumann entropy, 
and in the limit 9—7-00, which yields the distribution of the maximal eigenvalue. 
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